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1.  Introduction 


The  research  conducted  under  N00014-95-T0021  investigated  the  influence  of 
evaporating  sea-spray  under  high  wind  conditions.  These  investigations  were  acconplished  using 
a  interactive  Eulerian-Lagrangian  model  developed  by  the  principle  investigator.  During  the  first 
year  of  the  project,  the  interactive  model  was  validated  using  data  collected  during  the  1988 
Couche  Limite  UnidimensioneUe  Stationnaire  (CLUSE)  program.  These  efforts  were  reported  in 
“Spray  Droplet  Modeling,  II:  An  Interactive  Eulerian-Lagrangian  model  of  Evaporating  Spray 
Droplets”  by  Edson  et  al.  in  the  Journal  of  Geophysical  Research.  This  work  has  demonstrated 
that  the  combined  model  accurately  simulates  the  turbulent  transport  of  evaporating  droplets.  In 
additions,  this  paper  advanced  the  state-of-the-art  in  droplet  research  by  demonstrating  that  the 
potential  for  substantial  modification  of  the  surface  energy  budget  exists  if  the  presence  of  waves 
acts  to  eject  the  droplets  higher  and/or  disperse  the  droplets  more  efficiently. 

This  demonstration  was  taken  several  steps  further  in  an  abstract  presented  at  the  12”' 
Symposium  on  Boundary  Layers  and  Turbulence  entitled  “Modeling  the  Role  of  Sea  Spray  on 
Air-Sea  Heat  and  Moisture  Exchange”  by  Edson  and  Andreas.  This  abstract  compared  the  results 
from  the  interactive  model  with  a  model  developed  by  Dr.  Edgar  Andreas  of  the  U.S.  Army’s 
Cold  Regions  Research  and  Engineering  Laboratory.  The  agreement  between  these  models  was 
excellent  and  both  predicted  that  the  contribution  of  the  sea  spray  on  the  latent  heat  flux  is  at  least 
10%  of  the  total  under  high  wind  speed  conditions. 

The  second  part  of  this  investigation  has  involved  modifications  to  the  Eulerian  portion  of 
the  code  to  include  a  wavy  lower  boundary.  The  validation  of  this  model  is  being  acconq)lished 
through  conq)arisons  with  an  open  ocean  data  set  collected  aboard  the  R/P  FLIP  during  the 
second  part  of  the  1995  Marine  Boundary  Layers  (MBL)  Main  Experiment.  The  preliminary 
analyses  of  these  data  have  generated  a  great  deal  of  interest  in  the  modeling  community, 
particularly  because  they  questions  a  number  of  earlier  simulations  that  resulted  in  a  wave 
boundary  layer  (WBL)  that  extended  only  a  few  meters  above  the  surface.  Our  observations 
indicate  that  the  wave’s  influence  often  extended  beyond  the  height  of  our  upper  most  sensors  at 
18  m  A  sanqjle  of  these  observations  and  the  some  of  the  analysis  techniques  we  are  developing 
to  explain  them  are  given  sections  4  and  5. 
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These  observations  have  led  the  PI  to  develop  new  boundary  conditions  and  closure 
schemes  to  simulate  the  flow  over  waves  as  realistically  as  possible.  These  schemes  are 
summarized  in  the  final  section,  which  describes  the  PTs  work  in  progress.  The  goal  of 
this  work  is  to  validate  these  schemes  with  the  FLIP  observations.  The  coupled  model 
with  them  to  examine  the  influence  of  waves  on  sea  spray  dispersion. 
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Abstract.  This  paper  describes  an  interactive  Eulerian-Lagrangian  model  of  the 
turbulent  transport  of  evaporating  droplets.  A  k-z  (where  k  is  turbulent  Idnetic  energy 
and  8  is  its  rate  of  dissipation)  turbulence  closure  model  is  used  to  accurately  simulate 
stable,  near-neutral,  and  unstable  boundary  layers  within  die  large  air-sea  interaction 
tunnel  at  the  Institut  de  Mficanique  Statistique  de  la  Turbulence  (IMST),  Luminy, 
France.  These  results  are  then  used  with  the  Lagrangian  model  described  in  part  1 
[Edson  and  Fairall,  1994].  The  coupled  model  is  shown  to  give  excellent  agreement 
with  droplet  dispersion  measurements  made  during  the  1988  Couche  Limite 
Unidimensionelle  Stationnaire  d’Embrums  (CLUSE,  a  French  acronym  that  translates 
to  one-dimensional  stationary  droplet  boundary  layer)  campaign.  Additionally,  this 
paper  describes  how  the  coupled  model  can  now  be  used  to  investigate  the  interaction 
between  the  evaporating  droplets  and  the  turbulent  fields  of  temperature  and  humidity. 
The  investigation  shows  that  although  the  influence  of  the  droplets  is  small  under  the 
conditions  simulated  at  IMST,  the  potential  for  substantial  modification  of  the  surface 
energy  budget  exists  for  high-wind  conditions  over  the  ocean. 


1.  Introduction 

This  paper  describes  an  interactive  Eulerian-Lagrangian 
model  of  the  turbulent  transport  of  evaporating  spray  droplets. 
Ihe  model  has  been  developed  to  address  some  of  the  limita¬ 
tions  described  by  Edson  and  Fairall  [1994]  (hereinafter 
referred  to  as  part  1),  and  to  allow  the  use  of  the  model  in  more 
complicated  flows.  The  model  development  involved  the 
integration  of  the  Lagrangian  model  described  in  part  1  with  an 
Eulerian  model  of  turbulent  flows  that  uses  prognostic  equa¬ 
tions  for  the  evolution  of  the  turbulent  kinetic  energy  k  and  its 
rate  of  dissipation  € .  The  integrated  code  has  been  christened 
Gwaihir,  and  we  shall  refer  to  the  model  as  such  in  the  follow¬ 
ing  discussion. 

The  initial  tests  of  the  k-e  model  are  conducted  through 
simulations  of  developing  boundary  layers  using  a  two-dimen¬ 
sional  version  of  the  code.  The  model  runs  are  initialized  and 
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compared  with  measurements  taken  within  the  large  air-watCT 
interaction  simulation  tunnel  at  the  Institut  de  M6canique 
Statistique  de  la  Turbulence  (IMST),  Luminy,  France,  during 
the  1988  Couche  Limite  Unidimensionelle  Stationnaire 
dEmbrums  (CLUSE,  a  French  acronym  that  translates  to  one¬ 
dimensional  stationary  droplet  boundary  layer)  campaign 
[Mestayer  et  al.,  1990].  These  simulations  have  provided  a 
means  to  test  the  various  droplet  dispersion  aspects  of  Gwaihir, 
as  well  as  the  performance  of  the  Eulerian  code  in  simulations 
of  the  marine  atmospheric  surface  layer. 

The  paper  describes  in  some  detail  both  the  physical  model 
and  the  numerical  procedure  used  in  our  approach.  It  also 
addresses  some  of  the  advantages  of  this  combined  (Eulerian 
plus  Lagrangian)  approach  over  separate  approaches  (Eulerian 
or  Lagrangian)  in  simulations  of  the  turbulent  transport  of 
heavy  particles.  It  then  concludes  with  the  results  from  the 
interactive  model  for  simulations  of  droplet  dispersion  in  both 
a  laboratory  and  marine  atmospheric  surface  layer. 

2.  Eulerian  k-e  Model 

The  Eulerian  code  used  in  this  simulation  is  derived  from  2ik-e 
model  developed  at  the  Laboratoire  de  Mdcanique  des  Fluids, 
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Ecole  Centrale  de  Nantes,  France,  to  simulate  flows  around 
urban  structures  \L6vi  Alvaris  et  aL^  1990;  Livi  Alvar  is  and 
Sini,  1992;  Lakehal  et  aL,  1996].  In  the  atmospheric  surface 
layer,  expressions  for  the  instantaneous  velocity  field  for 
incompressible  fluid  flow  can  be  written 

dU, 


dU,  dU. 

— i  +U- — - 

dt  ^  dXj 


1  dP 

—-^-Si- 

Pa 


0  -0" 

V  V 

0: 


where  Einstein's  summation  notation  is  used  and  the 
Boussinesq  approximation  has  been  applied,  v  is  the  kinematic 
viscosity;  0^  is  the  ambient  virtual  potential  temperature; 
gj  =  (0,0,-^)  where  g  is  thegravitationd  acceleration;  and 
and  0y  and  virtual  potential  tanperature  of  the  reference  state 
of  the  fluid,  respectively  [Landahl  and  Mollo-Christensen, 
1986].  In  (2)  the  pressure  field  P  represents  the  departure 
from  the  reference  pressure  field  in  hydrostatic  balance.  In 
developing  equations  designed  to  study  flows  where  the  mean 
departure  from  hydrostatic  equilibrium  can  be  nonzero  (e.g., 
around  a  building),  Sim  [1986]  and  Sini  and  Dekeyser  [1987] 
decomposed  this  departure  from  hydrostatic  equilibrium  into 
mean  and  fluctuating  parts.  The  Reynolds  averaged  equations 
for  the  mean  variables  are  then  given  by 

du. 

^=0  (3) 

OXj 


2,1,  Closure 

The  Reynolds-stress  tensor  is  modeled  using  the  Boussinesq 
eddy  diffiisivity  concept 

_  dW:  dH]  2. 


where  is  the  Kronecker  delta  tensor,  \j  is  the  eddy  diffus- 
ivity,  and  k  is  the  turbulent  kinetic  energy  (TKE)  defined  as 

k  =  (8) 

Similarly,  the  scalar  fluxes  are  modeled  using 

(9) 


'  ’ax 


where  Kq  and  are  the  exchange  coefficients  for  potential 
temperature  and  specific  humidity,  respectively.  The  exchange 
coefficients  are  parametrized  as 

=  (11) 
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Kq  =  PrjVj. 
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where  is  a  model  coefficient  and  the  Prandtl  and  Schmidt 
munbers  for  turbulent  diffusion  are  assigned  the  same  value, 
Prj.=ScT  =0.95  [Hogstrom,  1988]. 

Closure  is  then  accomplished  through  prognostic  equations 
for  both  the  turbulent  kinetic  enwgy  and  its  rate  of  dissipation 

dk  _ 

Tt  ^dx. 

—  r  1  (14) 
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where  the  ovrbar  represents  an  ensemble  average;  lower  case 
letters  denote  the  turbulent  fluctuations;  the  total  specific 
humidity  Q  has  since  been  added  to  allow  for  the  inclusion  of 
moisture  in  the  model  equations;  is  the  specific  heat  at 
constant  pressure;  and  5^^  and  5^  represent  source  terms  for 
SOTsible  heat  and  moisture,  respectively.  The  source/sink  terms 
are  discussed  in  detail  below. 
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where  €  is  the  rate  of  irreversible  dissipation  of  kinetic  energy 
into  thermal  energy  and  o^,  and  are  model 

coefiBcients.  By  parameterizing  the  transfer  coefficients  in  this 
feshion,  we  believe  that  the  diffusive  properties  of  the  flow  are 
more  representative  of  the  intensity  k  and  persistence  kje  of  the 
turbulence  than  could  be  obtained  through  the  normal  applica¬ 
tion  of  first-order  closure  (i.e.,  /^-theory).  The  two-equation 
approach  is  also  less  complicated  and  less  expensive  to  run 
than  hi^er-order  models,  while  still  maintaining  many  of  their 
advantages  over  those  of  first  order  or  those  that  assume  a 
balance  between  production  and  dissipation  of  TKE. 

2.2.  Model  Coefficients 

The  numerical  coefficients  are  chosen  such  that  they  are 
representative  of  the  laboratory  boundary  layer  (LBL)  at  IMST. 
The  values  of  the  coefficients  andC^  are  based  on 

semiempirical  relationships  [e.g.,  Avva  et  al.,  1988; 
Duynkerke,  1988].  These  relationships  have  been  shown  to  be 
in  good  agreement  with  observation  made  within  the  air-sea 
interaction  tunnel  at  IMST  [Mestayer,  1980].  In  particular,  the 
near-surface  dissipation  rate  is  well  defined  in  neutral  condi¬ 
tions  by 


KZ 


where  is  the  velocity  scaling  parameter  (or  more  commonly 
the  friction  velocity)  and  the  von  Karman's  constant,  k,  is 
assigned  a  value  of  0.4  [Hogstroniy  1988]. 

We  use  the  measurements  of  Mestayer  [1980]  to  relate  the 
friction  velocity  to  vertical  velocity  variance  as  a^=  1.69w, . 
The  constant  of  proportionality  used  in  (1 1)  to  define  the 
eddy  diffusivity,  is  assigned  the  value  0.09.  This  value  has 
been  successfully  used  to  simulate  a  number  of  laboratory 
flows  [e.g.,  Launder  and  Spalding,  1974].  The  value  of 
is  found  from  measurements  taken  in  homogeneous  grid 
turbulence  where  the  diffusion  and  production  of  TKE  are 
negligible.  This  leads  to  a  situation  where  there  exists  a 
balance  between  the  advection  and  dissipation  of  TKE  such 
that  becomes  the  only  constant  of  significance  in  (14)  and 
(15).  Values  of  determined  from  these  measurements  are 
found  to  lie  within  the  range  of  1.8  to  2.0  [Avva  et  al.,  1988]; 
we  have  chosen  the  commonly  used  value  of  1.92  for  the 
present  model.  In  highly  stratified  flows,  where  the  Richardson 
number  has  reached  its  critical  value  Ri^,  the  transport  pro¬ 
cesses  again  become  negligible  and  (14)  and  (15)  can  be 
combined  to  give  [Duynkefke,  1988] 

c,i=Q(i-i?g  (17) 


Additionally,  near  the  surface,  where  we  expect  a  balance 
between  mechanical  production  and  dissipation  of  TKE,  a 
logarithmic  velocity  profile,  and  negligible  transport,  (15) 
reduces  to 


r  - 

^€l  ^€2 


1/2 


(18) 


Table  1,  Numerical  Constants  Used  in  the  Model 
Simulations  _ 


Constant 

Value 

C 

0.09 

P 

1.44 

1.92 

1.00 

o 

1.11 

€ 

K 

0.40 

PTj 

0.95 

SCrp 

0.95 

See  section  2.1  for  variable  definitions. 


Equations  (17)  and  (18)  can  then  be  combined  to  obtain 


CI^C^RK 


(19) 


In  the  present  model  we  have  assigned  the  values  1.0  and  0.25 
to  and  Ri^,  respectively.  Using  these  values  and  the  above 
expressions  we  obtain  the  values  of  the  coefficiaits  listed  in 
Tablet. 


23.  Numerical  Formulation 

The  numerical  formulation  is  an  adaptation  of  the  code 
Chensi  developed  by  Livi  Alvaros  [1992]  and  Levi  Alvar  is 
and  Sini  [1992]  for  an  inhomogeneous,  thiee-dimamonal  (3- 
D)  grid.  The  solution  of  the  system  of  equations  is  found  using 
the  Marker-and-Call  computing  method  as  presented  by  Hirt 
and  Harlow  [1967].  The  numerical  model  utilizes  a  staggered 
grid  configuration  that  defines  the  velocity  componoits  at  the 
cell  faces  and  the  scalar  variables  at  the  cell  centers  as  shown 
in  Rgure  1.  Variable  grid  spacing  is  used  to  allow  for  smaller 
grid  increments  in  regions  where  strong  gradients  are  expected, 
in  an  effort  to  reduce  the  numerical  noise  in  the  algorithm.  The 
dmvatives  are  determined  with  a  second-order,  finite  volume 
scheme  that  takes  into  account  the  variable  grid  ^cing  as 
desCTibed  by  Uvi  Alvaris  [1992]. 

The  numoical  method  is  explicit  in  time  and  uses  an  upwind 
weighted  difference  scheme  for  the  advection  terns  and  a 
centered  difference  scheme  for  the  diffusion  tenns.  The 
continuity  equation  is  satisfied  for  the  mean  velodties  at  each 
time  step  through  use  of  the  artificial  compressibility  method 
described  by  Chorin  [1967].  The  system  of  equations  is  then 
marched  forward  in  time  until  the  desired  level  of  convergence 
is  reached  such  that  the  Eulerian  variables  represoit  the  steady 
state  solution.  In  the  current  model  we  assume  that  a  steady 
state  has  been  reached  when  changes  in  model  parameters 
between  successive  time  steps  approach  the  computer’s 


5 


EDSON  ET  AL.:  EULERIAN-LAGRANGIAN  MODEL  OF  EVAPORATING  SPRAY 


T 

Azi 

1 


T 

ZZi 

1  i 


T 

Azzi 

1 


Zl  ZZi 


pi 

K 

B 

m 

II 

B 

B 

n 

H 

B 

a 

1 

1 

1 

i 

Figure  1.  Tlie  staggered  grid  configuration  used  in  the 
Eulerian  model 


numerical  precision  (currently,  a  Sun  Sparcstation  2  using 
double  precision). 

In  the  simulations  that  follow  the  height  of  the  lowest  grid 
point  is  set  to  1  cm.  Therefore  it  is  safe  to  assume  that  molecu¬ 
lar  effects  can  be  neglected  since  this  height  is  at  least  an  order 
of  magnitude  larger  than  the  roughness  lengdi  The  2-D 
domain  of  the  model  simulations  measures  50  m  by  0.85  m, 
which  is  rou^y  the  length  and  height,  respectively,  of  the 
turbulent  boundary  layer  at  IMST,  The  determination  of  the 
boundary  conditions  using  the  configuration  shown  in  Figure  1 
is  discussed  in  some  detail  below. 


where  the  subscripts  o  and  2  represent  the  values  at  z^  and  at 
the  second  grid  point,  respectively.  Equation  (20)  reduces  to 
the  classic  logarithmic  profile  when  a  constant  stress  layer 
prevails  near  the  surface;  that  is,  we  assume  that 

=  (21) 


The  roughness  length  at  the  end  of  the  tunnel  is  estimated 
from 

C 

V  /  ,2 

+  ^  (22) 


^^1/4,  1/2 


where  the  denotes  the  value  at  the  outflow  section  of  the 
model  domain.  The  first  term  on  the  right-hand  side  of  (22) 
gives  the  roughness  length  for  a  smooth  surface  with  C^=  9, 
while  the  second  term  is  based  on  Chamock's  formula 
[Charnock,  1955]  using  =  0.017  [Garratt,  1977].  The 
value  of  z^  at  the  lowest  ^d  point  in  the  inflow  section  is 
assumed  to  equal  that  for  a  smooth  surface.  Using  these  two 
values,  we  then  assume  that  the  surface  roughness  increases  in 
a  linear  fashion  with  fetch.  This  assumption  is  consistent  with 
actual  measurements  made  within  the  tunnel. 

The  vertical  velocity  at  the  lowest  grid  points  is  then  deter¬ 
mined  by  requiring  that  the  continuity  equation  is  satisfied  in 
the  lowest  cells.  First,  the  average  horizontal  velocity  through 
each  grid  face  Uj  is  found  by  integrating  the  log  profile  from  z^ 
to  Zj .  The  vertical  velocity  is  then  found  using  (3),  with  W^= 
0,  such  that 


Ax 


(23) 


where  Ajc  is  the  width  of  each  cell. 

A  condition  of  zero  difi&isive  flux  of  TKE  is  used  between  the 
first  two  grid  points 


2.4.  The  Boundary  Conditions 

In  the  k-e  code  we  define  the  grid  as  shown  in  Figure  1, 
where  the  parameters  in  the  first  grid  are  not  defined  at  its 
center.  This  approach  saves  computer  time  and  improves  the 
determination  of  the  derivatives  at  the  lowest  grid  points.  It 
also  allows  us  to  move  the  height  of  the  first  grid  point  to  a 
re^on  well  above  the  viscous  sublayer  such  that  we  can  safely 
ignore  molecular  effects.  However,  these  noncentered  grid 
points  found  in  the  first  cells  require  special  consideration.  The 
velocity  component  parallel  to  the  surface  (in  these  simulations 
the  horizontal  velocity)  is  found  using  the  wall  law  given  by 
Launder  and  Spalding  [1974]  to  estimate  the  surface  stress 


The  dissipation  rate  of  TKE  is  found  by  assuming  that  the 
mechanical  production  of  TKE  is  equal  to  dissipation  in  the 
near-wall  region.  This  leads  to  a  relationship  between  the  TKE 
and  its  dissipation  rate  at  the  lowest  grid  point  given  by 


The  horizontal  velocity  profile  at  the  upstream  boimdary 
(hereafter  referred  to  as  the  entrance)  is  fixed  using  a  profile 
based  on  the  measurements  made  by  Selva  [1979]  in  the  IMST 
tunnel  using  a  Pitot  tube  and  highly  accurate  manometer. 
Neumann  conditions  (i.e.,  zero  gradient)  are  used  for  all 
variables  at  the  upper  and  downstream  boundaries,  except  for 
the  vertical  velocity  at  the  upper  boundary,  which  is  set  equal 
to  zero. 
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The  initial  values  of  k  are  found  by  inversion  of  (21),  while 
initial  values  of  e  are  assigned  using  (16).  The  value  of  w, 
used  in  both  relationships  is  assigned  an  initial  value  equal  to 
two  fifths  the  value  of  measured  at  a  distance  (fetch)  of  30 
m  from  the  tunnel  entrance.  The  inflow  values  of  kmde  are 
then  allowed  to  adjust  through  the  use  of  Neumann  conditions. 
The  temperature  and  specific  humidity  profiles  at  the  tunnel 
entrance  are  fixed  using  the  value  measured  at  0.75  m,  except 
for  at  the  lowest  grid  point  that  is  given  the  surface  value.  This 
is  consistent  with  both  measurements  and  the  action  of  the  heat 
exchangers  in  the  tunnel,  which  act  to  hold  the  air  temperature 
and  dew  point  constant,  mix  the  air  thoroughly,  and  break  down 
any  undesirable  eddies.  At  the  lower  boundary  the  temperature 
is  givCTi  the  value  of  the  water  surface  temperature  while 
the  specific  humidity  is  assigned  its  saturation  value  at 

2.5.  Eulerian  Model  Results 

The  Eulerian  model's  velocity  results  are  in  good  agreement 
with  the  measurements  made  using  a  Pitot  tube  during  the 
CLQSE  campaign  as  shown  in  Figure  2.  The  various  curves  in 
Figure  2  depict  the  evolution  of  the  wind  profile  at  various 
fetches.  The  velocity  profile  measured  at  30  m  gives  excellent 
agreanent  with  the  model  profile  at  a  fetch  of  29  m.  This  result 
indicates  that  the  velocity  evolution  is  accurately  modeled, 
especially  since  there  is  some  uncertainty  as  to  where  we  define 
the  fetch  to  be  equal  to  0.  The  small  discrepancy  between  the 
model  and  measurements  at  the  top  of  the  boundary  layer  is  the 
result  of  the  confluence  of  the  tunnel's  two  boundary  layers. 

In  general,  the  vertical  structure  of  temperature  and  humidity 
along  the  length  of  the  tunnel  cannot  be  adequately  described 


Figure  2.  A  comparison  of  the  velocity  measurements  made 
at  a  fetch  of  30  m  in  the  Institut  de  M6chanique  Statistique  de 
la  Turbulence  (IMST)  tunnel  with  the  Eulerian  model  results. 
The  group  of  curves  d^icts  the  model-derived  wind  profiles  at 
the  fetches  indicated.  The  error  bars  denote  the  standard 
deviation  of  the  measured  velocity  variance. 


Temperature  C  C) 


Figure  3.  A  comparison  of  the  temperature  measurements 
made  at  a  fetch  of  30  m  in  the  IMST  tunnel  with  the  Eulerian 
model  results.  The  three  sets  of  data  are  representative  of  the 
high-,  medium-,  and  low-  humidity  (95%,  77%,  and  55%, 
respectively)  runs  made  during  CLUSE.  The  nominal  wind 
speed  in  all  three  cases  is  approximately  7.5  m  s'^  The  group 
of  curves  depicts  the  model-derived  temperature  profiles  at  the 
fetches  indicated.  The  standard  deviation  of  the  measured 
temperature  variance  is  given  at  each  level  for  the  medium 
humidity  run. 


-1 

Specific  Humidity  (gKg  ) 


Figure  4.  A  comparison  of  the  specific  humidity  measure¬ 
ments  made  at  a  fetch  of  30  m  in  the  IMST  tunnel  with  the 
Eulerian  model  results.  The  conditions  are  the  same  as  in 
Figure  3.  The  standard  deviation  of  the  measured  specific 
humidity  variance  is  giv^  at  each  level  for  the  medium 
humidity  run. 
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Figure  5,  A  comparison  of  the  momentum  flux  measured  at 
the  indicated  wind  speeds  at  a  fetch  of  30  m  in  the  IMST  tunnel 
with  the  Eulerian  model  results.  The  group  of  curves  depicts 
the  growth  of  the  simulated  boundary  layer  through  flux 
profiles  determined  at  27, 29, 31,  and  33  m  (from  right  to  left, 
looking  at  the  top  of  the  curves)  for  a  wind  speed  of  7.5  m  s‘\ 


by  single  profiles  as  used  in  part  1  (even  without  the  droplets). 
Owing  to  the  mixing  of  the  air  by  the  turboprop  and  heat 
exchangers,  the  air  enters  the  test  section  with  approximately 
constant  values  of  specific  humidity  and  temperature.  As  the 
air  moves  along  the  length  of  the  tunnel,  it  interacts  with  the 
water  surface  through  molecular  and  eddy  exchange.  This 
interaction  is  shown  in  Figures  3  and  4,  which  indicate  good 
agreement  between  measured  and  modeled  results  at  a  fetch  of 
30  m.  These  measurements  were  made  using  a  thermocouple 
psjchrometer  system  provided  by  the  University  of  Washington 
using  chromel  constantan  thermocouples.  These  devices  are 
believed  to  be  accurate  to  ±  0.2®C;  however,  cold  spikes 
caused  by  droplets  impacting  the  sensors  are  expected  to 
increase  the  uncertainty  in  these  measurements  [Katsaros  et 
al,  1994]. 

The  result  of  primary  interest  for  droplet  dispersion  is  shown 
in  Figure  5,  which  displays  the  -uw  component  of  the 
Reynolds-stress  tensor  computed  from 

=  (26) 

oz 

with  carefully  conducted  measurements  made  at  a  fetch  of  30 
m  over  a  range  of  wind  speeds  by  Giovanangeli  and  LeCalve 
[1990]  at  IMST.  The  group  of  curves  representing  the  model 
results  at  several  fetches  shows  that  the  model  gives  good 
agreement  with  the  measurements  made  at  the  same  wind 
speed  at  a  fetch  of  30  m.  These  curves  are  representative  of  a 
developing  LBL  and  indicate  that  the  growing  surface  layer 
reaches  a  height  of  approximately  0.2  m  at  30  m.  The  implica¬ 
tions  of  this  developing  surface  layer  on  droplet  dispersion  is 
addressed  in  section  3. 


3.  Droplet  Dispersion  Modeling 

Several  Eulerian  approaches  have  been  successfully  applied 
to  the  study  of  turbulent  diflusion  of  discrete  particles  in 
boundary  layer  flows.  These  include  the  studies  of  Ling  and 
Kao  [1976],  Ling  et  al  [1978, 1980],  Burk  [1984],  Mostafa 
and  Mongia  [1987],  Stramska  [1987],  and  Rouault  et  ai 
[1991],  among  others.  However,  the  restraints  placed  on  these 
models  often  restrict  their  application  to  either  very  small, 
highly  concentrated  particles  or  simple  flow  geometries  with  a 
homogeneous  source  of  particles  that  can  then  be  studied  using 
budget  equations. 

Note  that  these  constraints  are  not  always  a  drawback.  For 
example,  in  highly  concentrated  flows,  where  the  particle- 
particle  interactions  are  not  negligible,  it  is  g^erally  much 
easier  to  include  these  effects  in  a  Eulerian  scheme  than  a 
Lagrangian  one.  Additionally,  the  use  of  budget  equations  can 
provide  a  means  of  studying  the  effects  of  particle  interaction 
with  the  scalar  fields  through  the  use  of  source/sink  relation¬ 
ships,  which  are  difficult  to  include  in  a  purely  Lagrangian 
scheme.  This  approach  has  been  successfully  employed  by 
Rouault  et  al  [1991],  who  studied  the  effect  of  droplet 
evaporation  on  the  scalar  fields  of  temperature  and  humidity. 

Tlie  simidation  of  the  movement  of  heavy  particles  involves 
parameterizations  to  account  for  the  effects  of  gravity  and 
inertia.  These  parameterizations  gaierally  require  the  tuning 
of  adjustable  constants  through  comparison  with  experimental 
data  [e.g.,  Rouault  et  aL,  1991].  Unfortunately,  these  calibra¬ 
tion  measurements  are  difficult  to  come  by,  several  notable 
exceptions  being  the  studies  of  Snyder  and  Lumley  [197 1]  and 
the  Humidity  Exchange  Over  the  Sea  (HEXOS)  in  a  Simula¬ 
tion  Tunnel  (HEXIST)  experiments  [e.g.,  Andreas  et  ai, 
1 995] .  There  is  also  a  question  as  to  the  universality  of  these 
constants  when  they  are  applied  to  more  complicated  flows. 
iWilly,  as  in  any  Eulerian  model,  discrete  sources  of  particles 
are  difficult  to  include  in  the  calculation  domain  and  generally 
require  subgrid-scale  parameterizations  that  are  often  difficult 
to  quantify. 

Therefore  we  feel  that  the  Lagrangian  approach  is  the  best 
alternative  if  one  is  concerned  with  the  dispersion  of  low- 
concentration,  heavy  particles  from  discrete  or  nonuniform 
sources  (e.g.,  from  spume  droplets  produced  at  the  wave 
crests)  and  in  fluid  flows  with  complicated  geometries  (e.g., 
over  waves  or  within  a  surf  zone).  Tb&  approach  also  allows 
for  much  more  flexibility  in  specifying  boundary  conditions. 
For  example,  a  Lagrangan  scheme  can  easily  include  boundary 
conditions  where  the  particles  either  rebound  or  stick,  depend¬ 
ing  on  the  type  of  surface  one  is  trying  to  simulate.  It  is  also 
especially  advantageous  in  situations  where  the  physicochemi¬ 
cal  characteristics  of  the  particles  change  rapidly  when  they 
experience  a  highly  inhomogeneous  environment.  In  this 
instance,  time  rate  of  change  equations,  which  may  depend 
upon  the  local  conditions,  as  well  as  on  the  particle’s  history, 
are  easily  included  in  the  model  as  long  as  these  equations  are 
known. 

The  Lagrangian  model  is  described  in  detail  in  part  1. 
Briefly,  the  model  relies  on  a  finite  difference  form  of  the 
Langevin  equation  for  the  fluctuating  components  of  the 
droplet's  velocity  due  to  turbulence.  When  the  mean  droplet 
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fall  velocity  is  added  to  this  fluctuating  component  w'’(0, 
the  instantaneous  vertical  velocity  is  given  by 


WP(t+At)  = 


‘  A 

1-^ 


WP(S) 


w 


(27) 


whore  is  the  droplet  integral  timescale,  Ar  is  the  time  step, 

is  the  ^uare  root  of  the  droplet  vertical  velodty  variance,  and  C(0 
is  a  random  number  drawn  from  a  Gaussian  distribution  with 
zero  mean  and  variance  At.  The  expressions  for  and 
are  then  derived  from  the  equation  of  motion  for  small,  heavy 
droplets  given  by 


=  (28) 

Dt  a  a 


where 


1  + 


(32) 


Hie  parametiM*  1/A  accounts  for  the  decorrelating  effect  of  the 
droplet  falling  out  of  a  fluid  parcel  where  the  fluctuating  fluid 
velocities  are  highly  correlated  over  small  time  steps.  In  past 
studies  this  parameter  has  been  tuned  to  fit  experimental  data 
by  multiplying  (32)  by  a  constant  ranging  between 

0.3  and  1  [e.g.,  Mostitfa  and  Mongia,  1987;  Wells  and  Stock, 
1983],  In  this  study  we  continue  to  use  the  formulation  given 
in  part  1;  that  is,  the  constant  is  left  equal  to  1.  We  address  the 
implication  of  this  choice  in  the  discussion  given  below. 

The  Lagrangian  integral  timescale  is  commonly  defined  in 
engineering  literature  [e.g,,  Mostcfa  and  Mongia,  1987]  using 
the  output  of  the  k-e  model  as  Xj^  =  ^k/e,  where  P  ranges 
between  0.15  and  0.6.  If  we  equate  this  expression  with  the 
expression  used  in  part  1,  we  obtain 


where  the  superscripts  /and  R  denote  the  fluid  and  relative 
velocities,  respectively,  and  the  total  derivative  denotes  motion 
following  the  droplet  The  parameter  K  represents  the  ratio  of 
the  Stokes  velocity  to  the  mean  relative  fall  velocity,  and  a  is 
the  response  time  for  droplets  that  ideally  obey  Stokes  law. 

The  approach  given  in  part  1  yields  the  following  expression 
for  the  vertical  velocity  variance: 


(1  +x) 


(29) 


The  vertical  velocity  variance  is  parameterized  in  terms  of  the 
kinetic  energy  by  combining  (21)  with  the  results  from 
Mestayer  [1980]  to  obtain 

gI=  1.69  (34) 


4.  Gwaihir 


where 

X-^  (30) 


The  parameter  %  is  the  ratio  of  the  droplet  response  time  to 
flie  T  agrangian  integral  time  scale  x^.  This  parameter  thereby 
determines  how  the  droplet  reacts  to  die  turbulent  motion  of  the 
surrounding  fluid,  e.g.,  as  the  droplet  encounters  smaller  eddies 
as  it  nears  the  surface  (smaller  x^),  the  iirfluence  of  the 
turbulence  on  the  droplet  motion  diminishes  because  the 
droplet  can  no  longer  react  to  these  smaller  eddies  due  to  its 
inertia.  The  height  at  which  this  begins  to  occur  is  determined 
by  the  droplet’s  size  as  reflected  by  its  response  time  a/K,  The 
equivalent  situation  occurs  as  the  turbulence  intensity  increases 
(again,  smaller  x^)  such  that  its  response  time  becomes  too 
large  to  permit  the  droplet  to  react  fast  enough  to  all  scales  of 
the  turbulence. 

The  integral  timescale  is  a  measure  of  the  persistence  of  the 
droplet’s  velocity  as  it  moves  through  the  fluid.  This  coefficient 
is  derived  by  detennining  the  droplet  autocorrelation  coeffident 
in  an  ^proach  similar  to  that  used  by  Meek  and  Jones  [1973], 
The  droplet  integral  timescale  is  then  determined  by  integrating 
the  autocorrelation  function  over  time,  which  yields 

<  =  ^(ln)  (31) 


With  these  parameters  we  have  all  the  necessary  expressions 
to  simulate  the  turbulent  dispersion  of  heavy  particles  once  the 
Eulerian  fields  are  available  from  the  k-e  code.  Therefore, 
z&er  a  steady  state  solution  of  a  particular  flow  configuration  is 
found,  the  velocity  and  scalar  values  are  passed  to  tiie 
Lagrangian  section  of  the  code  in  order  to  detonine  the 
parametes  in  the  above  expressions.  These  parameterizations 
are  then  used  with  (27)  to  compute  the  trajectories  of  evaporat¬ 
ing  drqilets.  The  amount  of  trajectories  begun  over  a  specified 
time  T  is  determined  by  a  usa'-defined  source  function.  As 
these  trajectories  are  computed,  the  Lagrangian  code  keeps 
track  of  the  droplet  concentration  ,  sensible  heating  rate  S^, 
and  water  vapor  production  in  each  cell.  Once  the 

Lagrangian  portion  of  the  code  has  ejected  all  of  the  droplets 
determined  by  the  source  function  over  the  time  T,  tiie 
Eul^ian  portion  of  the  code  is  then  renm  with  the  new  nonzero 
values  of  qp,  and  Sm-  This  process  is  repeated  until  the 
k-e  code,  upon  completion  of  the  last  Lagrangian  run,  reaches 
a  steady  state  over  an  interval  less  than  T. 

The  source  of  droplets  is  determined  from  the  surface  source 
fiinction,  which  gives  the  number  of  droplets  per  unit  time,  per 
unit  area  produced  at  the  surface.  Since  the  droplet  measuring 
devices  were  situated  slightly  downwind  of  the  30-m-long 
bubbler  array  used  during  the  CLUSE  campaign  [Mestayer  et 
aL,  1990],  the  Gwaihir  source  function  was  determined  by 
matching  the  model  output  with  droplet  measurements  made  at 
a  height  of  0.2  m,  0.5  m  downwind  of  the  simulated  source. 
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The  total  number  of  droplets  over  a  radius  interval  between 
r-Ar/2  and  r+Ar/2  ejected  from  a  prescribed  area  (where  Ar 
is  the  radius  increment)  is  given  by 

N(r)  =  Nj,(r)T:^rAx  Ay  (35) 


where  Nj(r)  is  the  surface  source  function,  Tis  the  total  time 
that  the  source  is  on,  Ax  is  the  length  of  the  area,  and  Ay  is  its 
width  (set  equal  to  unity  in  this  2-D  simulation).  Since  Gwaihir 
uses  a  horizontal  ^d  spacing  of  Ax  =  0,5  m,  the  source  is 
divided  up  into  60  sections. 

The  droplets  produced  in  each  section  are  then  released  at 
their  ejection  height.  This  height  corresponds  to  the  height  to 
which  a  jet  droplet  is  ejected  by  a  bubble  bursting.  The  height 
and  number  of  droplets  ejected  are  a  ftinction  of  bubble  size, 
and  we  have  used  the  data  given  by  Blanchard  and  Woodcock 
[1957]  to  determine  these  parameters  as  described  in  part  1. 

As  the  released  droplet  moves  through  the  model  domain,  the 
values  needed  in  the  Lagr^gian  model  parameterizations  are 
determined  from  a  linear  average  of  the  four  nearest  Eulerian 
grid  points  in  order  to  find  C/(X,Z) ,  W(X,Z),  and  $(X,Z)  at 
X,Z.  Owing  to  the  staggered  grid,  this  involves  keeping  track 
of  four  different  indices  corresponding  to  x,xx,  2,  and  zz.  If 
the  droplet  falls  below  zzj ,  the  values  of  the  variables  are  found 
using  the  same  wall  functions  employed  in  the  Eulerian  model. 
If  the  droplet  falls  below  z^,  it  is  assumed  to  stick  to  the 
surface,  in  which  case  another  trajectory  is  begun.  A  new 
droplet  trajectory  is  also  begun  if  a  droplet  is  carried  out  of  the 
calculation  domain. 

The  concentration  of  droplets  over  a  specified  size  range  is 
calculated  by  accumulating  the  time  these  droplets  spend  in 
each  cell.  The  accumulated  time,  divided  by  the  total  time  that 
the  source  is  on,  divided  by  the  volume  of  the  cell  is  a  measure 
of  the  droplet  concentration.  The  droplet  volume  concentration 
(droplet  volume,  per  unit  volume,  per  radius  increment)  in  each 
cell  is  then  computed  from 


dy(,r)_4  y 

dr  3VTAr^ 


(36) 


where  the  summation  is  for  each  cell  over  the  advection  time 
T,  and  V = AxAyAz  is  the  volume  of  each  cell. 

The  source  functions  of  p^c^Q  and  (2  are  the  means  by 
which  the  evaporating  ^oplets  interact  with  the  fields  of 
sensible  heat  and  moisture.  These  functions  can  act  as  either 
a  source  or  sink  of  sensible  h^t  and  moisture  depending  on  the 
ambient  conditions.  Therefore  Q  in  (6)  is  actually  the  total 
specific  humidity 

Qixa)  =  +  qjpca)  (37) 


In  breaking  down  the  specific  humidity  in  this  fashion  the 
source/sink  functions  become  solely  a  function  of  droplet 
evaporation/condensation.  The  interaction  between  the 
evj^orating  droplets  and  the  scalar  fields  is  then  simulated  by 
releasing  all  the  droplets  produced  during  a  given  period  of 
time  and  accumulating  the  sensible  heat  and  moisture  they 
consume  or  release  in  each  cell.  The  functions  used  in  this 
model  are  given  as 

(38) 


^  dt 


(39) 


where  the  summation  is  over  all  droplets  in  each  cell,  is  the 
density  of  the  droplet,  is  the  ventilation  coefficient  for  heat 
diffusion  [Beard  and  Pruppacher,  1971],  is  the  thermal 
conductivity  of  air,  and  and  are  the  air  and  droplet 
surface  temperatures,  respectively. 

The  ev^oradon  rate,  dridt,  is  determined  using  the  function 
given  by  Beard  and  Pruppacher  [1971]  using  modification 
proposed  by  Andreas  [1989]  for  curvature  effects,  while  the 
droplet  surface  temperature  is  governed  by  the  functions 
described  in  part  1.  The  droplet  surface  temperature  is 
assigned  the  value  of  the  wato*  surface  at  the  time  of  its  release. 
Initialized  in  this  way,  Gwaihir  explicitly  models  the  initial 
transfer  of  sensible  heat  from  droplet  to  air  when  the  sea 
surface  temperature  exceeds  the  air  temperature  (unstable 
conditions)  as  described  by  Andreas  et  al  [1995],  as  well  as 
the  eventual  cooling  (heating)  of  the  atmosphere  due  to 
evaporation  (condensation). 

The  virtual  ten^rature  r^uired  in  (4),  (14),  and  (15)  is  then 
determined  by  considering  all  three  sources  of  moisture  [Stull, 
1988] 

=  (40) 

0(x^)  [1  +  0.61  q^(x^)  -  q^(x^)  -  qpix^W  ’ 


where  is  the  contribution  to  the  total  specific  humidity  from 
the  droplets;  that  is,  it  is  found  by  integrating  the  droplet 
volume  concentration  in  each  cell.  In  practice,  is  equal  to 
zero  unless  Q  is  greater  than  its  saturation  v^ue.  If  siq)er- 
saturation  occurs,  then  ^/x,z)  is  equal  to  q^{x^)  and  ^^(x,z) 
is  equal  to  Q{x^)  -  q^j^x^) . 

These  additional  sources  of  moisture  are  included  in  the 
buoyancy  flux  using  the  approach  given  by  Stull  [1988] 


where  q^  is  the  specific  humidity  due  to  water  vapor  and  is 
the  contribution  due  to  liquid  water  other  than  the  spray 
droplets  (e.g.,  fog).  H^e  we  are  assuming  that  the  fog  droplets 
are  small  enough  that  they  can  be  properly  modeled  by  (6);  that 
is,  the  fog  droplets  are  small  enough  to  be  treated  as  passive 
scalars.  On  the  other  hand,  the  spray  droplets  are  too  large  for 
such  treatment  and  are  instead  modeled  separately  by  the 
Lagrangian  code. 


•  =  iv0  [1  +  0.61^^  -qj^-  q^ 

+  Q^.6\wq-wq^- 

where  we  are  ignoring  the  triple  products  found  in  its  full 
derivation.  The  fluxes  are  defined  by  breaking  down  the 
moisture  terms  as  described  above  and  using 
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-wC’  =  Pr^j^  (42) 

02 

where  and  we  assume,  following  Rouault  et  ql 

[1991],  that  the  turbulent  Prandtl  and  droplet  Schmidt  numbers 
are  equal. 

5.  Gwaihir  Results 

In  this  s^tion  we  concentrate  on  how  we  have  addressed 
some  of  the  earlier  model’s  deficiencies  by  combining  it  with 
the  Eulerian  code.  Specifically,  we  examine  the  improvement 
between  the  model  and  measurement  comparisons  due  to  a 
more  realistic  simulation  of  the  turbulent  field  through  which 
the  droplets  ^e  dispersing.  Second,  we  demonstrate  the 
in5)roveinents  in  our  simulations  due  to  the  ability  of  Gwaihir 
to  pmnit  interaction  between  the  droplets  and  the  scalar  fields. 
Finally,  we  examine  the  extent  to  which  the  droplets  modify  the 
scalar  fields  and  examine  how  this  effect  modifies  the  fluxes  of 
sensible  and  latent  heat, 

5.1.  Influence  of  the  Turbulent  Flow  Field  on  Droplet 
Dispersion 

In  section  2.5  we  demonstrated  that  the  model  is  capable 
of  accurately  simiilating  a  2-D  developing  surface  layer.  We 
now  examine  how  the  improved  simulation  of  the  turbulent 
flow  field  affects  the  output  of  the  Lagrangian  component.  In 
Figure  6,  we  begin  our  comparison  between  Gwaihir  and  the 
droplet  measurements  tak^  during  the  Q.USE  campaign.  It 
shows  excellent  agreement  between  Gwaihir  (at  a  fetch 
corresponding  to  30.5  m)  and  measurements  for  dl  sizes  and 


Figure  Droplet  volume  spectra  as  a  function  of  droplet 
radius.  The  solid  lines  represent  the  model  results,  while  the 
symbols  represent  data  measured  at  the  indicated  heights.  The 
error  bars  are  computed  at  every  fifth  data  point  from  dV/dr  x 
(1  ±  ll\/N)y  where  N  is  the  number  of  droplets  counted  over 
the  size  interval.  The  conditions  during  the  run  were  a  relative 
humidity  of  95%  and  a  nominal  wind  speed  of  7.5  m  s'\ 


Figure  7.  Profiles  of  total  liquid  water  content  q^^  ,  calculated 
by  integrating  over  the  measured  and  modeled  droplet  spectra 
at  each  height.  The  symbols  represent  data  measured  at  the 
indicated  heights.  The  solid  line  represents  Gwaihir  results 
taken  directly  over  the  simulated  source  of  droplets,  while  the 
dashed  line  shows  the  profile  1  m  beyond  the  end  of  this 
source.  The  dash-dotted  line  reproduces  the  one-dimensional 
model  results  given  by  Edson  and  Fairall  [1994].  The 
conditions  during  the  run  were  a  relative  humidity  of  95%  and 
a  nominal  wind  speed  of  7.5  m  s’*. 


heights,  except  for  the  largest  droplets  at  63  cm.  This  is 
perhaps  a  result  of  the  merging  of  two  boundary  layers  in  the 
actual  tunnel,  which  we  have  not  tried  to  model  in  this  simula¬ 
tion. 

Upon  close  inspection  of  Figure  6  we  see  that  Gwaihir  is 
slightly  underestimating  the  vertical  dispersion  of  the  larger 
droplets.  The  agreement  between  the  model  and  measurements 
at  the  largest  sizes  could  be  improved  by  including  a  fractional 
constant  in  (32)  as  described  above.  Howevo*,  its  inclusion 
would  worsen  the  agreement  at  smaller  sizes.  A  more  likely 
e)q)lanation  for  the  disagreement  is  that  we  are  simply  pushing 
Gwaihir  too  far;  that  is,  we  are  using  a  1-D  Lagrangian  model 
in  a  2-D  flow.  A  better  alternative  is  to  make  the  Lagrangian 
portion  of  the  code  2-D  by  including  a  correlated  horizontal 
component  of  the  droplet  velocity  as  did  Ley  and  Thomson 
[1983],  Additionally,  the  assumptions  that  lead  to  its  develop¬ 
ment  suggests  that  it  works  best  in  homogeneous  turbulence, 
which  is  clearly  not  the  case  in  the  LBL.  In  fact,  we  expect 
Gwaihir  to  work  better  in  simul^ons  of  a  marine  surface  layer, 
where  these  conditions  are  better  realized. 

Even  with  these  minor  shortcomings.  Figures  7  and  8  indicate 
that  we  have  improved  the  model  simulations  (compared  with 
the  results  given  in  part  1)  by  simply  impro-ving  the  simulation 
of  the  turbulent  flow  field,  rather  than  trying  to  tune  the 
Lagrangian  model  parameters.  This  is  evident  in  the  profiles 
of  q^  given  in  Figure  7,  which  is  computed  by  approximating 


11 


EDSON  ET  AL.;  EULERIAN-LAGRANGIAN  MODEL  OF  EVAPORATING  SPRAY 


D 


Figure  8.  Profiles  of  total  liquid  water  content  q^,  for  a 
relative  humidity  of  55%  and  a  nominal  wind  speed  of 
7.5  m  s'*.  The  lines  are  as  denoted  in  Figure  7. 


the  area  under  the  curves  shown  in  Hgures  6  and  9  in  part  1- 
using 


Pa  \  i 


where  N  is  the  total  number  of  radius  increment  bins  and  the 
density  units  are  chosen  to  give  grams  per  kilogram.  Figure  7 
shows  that  either  version  of  the  model  gives  excellent  agree¬ 
ment  with  the  measurements  within  the  constant  flux  layer. 
Recall  that  the  constant  flux  is  assumed  to  exist  throughout  the 
1-D  boundary  la>er  simulated  in  part  1 .  This  is  the  reason  why 
Gwmhir  gives  much  better  agreement  as  we  near  the  top  of  the 
boundary  layer,  where  the  momentum  flux  tends  toward  zero, 
which  causes  a  drastic  reduction  in  the  vertical  dispersion  of 
the  droplets. 

Figure  8  shows  that  Gwaihir  does  a  better  job  of  handling 
droplet  evaporation  than  the  l-D  model,  which  tended  to 
overevaporate  the  droplets.  This  is  due  to  the  inclusion  of 
droplet  feedback  mechanisms  (i.e.,  the  moistening  of  the  near 
surface  which  lessens  the  amount  of  evaporation),  as  well  as 
the  general  improvement  of  the  mean  profile  simulations  of 
temperature  and  moisture.  In  fact,  this  latter  effect  is  most 
likely  the  cause  for  the  improvement,  for  reasons  given  in 
section  5.2, 


5.2.  Influence  of  Droplets  on  the  Scalar  Fields 

The  first  step  in  the  simulation  is  to  compute  the  turbulent 
fields  of  velocity,  temperature,  and  humidity  in  the  absence  of 
droplets.  Therefore  the  k-6  model  gives  us  an  easy  way  to 
detamine  the  effect  of  droplet  evaporation  on  the  scalar  fields 
by  examining  the  diffaence  in  the  temperature  and  humidity 
profiles  modeled  with  and  without  droplets.  In  Figure  9  we 
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F^ure  9.  Profiles  of  the  change  in  the  model-drived  tempera¬ 
ture  profiles  due  to  the  presence  of  droplets.  The  nominal  wind 
speed  was  7.5  m  s''  for  all  runs,  while  the  relative  humidity  is 
as  indicated  on  the  profiles  measured  at  30  m  (55%,  77%,  and 
95%). 


depict  the  difference  in  temperature  profiles  due  to  the  presence 
of  droplets.  The  high  huniidity  run  in  Figure  9  illiistrates  a 
case  where  the  release  of  sensible  heat  firom  the  droplets  causes 
the  near-surface  air  temperature  to  actually  increase.  This  can 
be  attributed  to  the  large,  positive  water-air  temperature 
difference  and  little  evaporative  cooling  at  high  humidity. 


Specific  Humidity  Change  CgKg  ) 


Figure  10.  Profiles  of  the  change  in  the  model-derived  specific 
humidity  profiles  due  to  the  presence  of  droplets.  The 
nominal  wind  speed  was  7.5  m  s''  for  all  nins,  while  the 
relative  humidity  is  as  indicated  on  the  profiles  measured  at  30 
m  (55%  and  77%). 
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The  lower-humidity  runs  shown  in  Figures  9  and  10  give  the 
expected  result  of  lower  near-surface  temperatures  and  higher 
specific  humidities  in  the  presence  of  droplets  due  to  their 
evaporation.  The  magnitude  of  this  change  is  similar  to  the 
results  of  Rouault  etaL  [1991],  with  a  maximum  of  approxi¬ 
mately  0.1  in  temperature  and  0.05  g  k  g‘'  in  specific 
humidity.  This  is  a  barely  measurable  change  (in  fact,  the 
precision  of  the  sensors  may  not  be  able  to  resolve  it)  and  is 
smaller  than  the  measurements  shown  in  Figure  16  in  part  1. 
As  stated  in  part  1,  we  believe  that  this  is  a  consequence  of 
droplets  wetting  the  "dry-bulb"  thermometer  resulting  in  an 
underestimation  of  the  dry-bulb  temperature  and  overestima¬ 
tion  of  the  specific  humidity.  For  this  reason  we  believe  that 
the  improvement  in  the  model’s  performance  is  most  likely  due 
to  the  inclusion  of  the  2-D  simulation  of  the  mean  temperature 
and  humidity  profiles,  rather  than  the  droplet  feedback  mecha¬ 
nism. 

5,3,  Influence  of  the  Droplets  on  the  Surface  Energy 
Budget 

On  the  basis  of  the  results  given  in  section,  5.2,  one  might 
conclude  that  we  should  not  expect  droplets  to  have  an  appre¬ 
ciable  effect  on  the  surface  energy  budget.  However,  there  are 
a  number  of  reasons  why  this  may  not  be  the  case  for  the 
marine  boundary  layer  (MBL).  For  example,  we  have  already 
shown  that  the  intensity  of  mv  in  the  LBL  falls  off  rapidly  with 
height  This  decrease  has  a  particularly  adverse  effect  on  large 
droplet  dispersion  as  evidenced  by  Figures  7  and  8,  where  the 
majority  of  the  potential  moisture  from  the  droplets  is  found 
below  the  highest  height  of  0.18  m,  indicating  that  the  largest 
droplets  (which  dominate  the  total  water  content)  are  rarely 
found  above  this  level.  This  effect  is  also  responsible  for  the 
reduction  of  the  droplet  concentration  profile,  as  one  moves  0.5 
m  downwind  of  their  source,  due  to  the  rapid  fallout  of  the 
largest  droplets. 

In  fact,  direct  measurements  of  droplet  profiles  made  from  a 
wave  follower  by  DeLeeuw  [1986, 1987]  have  shown  that  the 
gradient  over  the  ocean  is  much  smaller  than  that  measured  in 
the  laboratory.  While  the  mechanisms  responsible  for  the 
enhanced  mixing  remain  a  hotly  debated  topic,  it  has  been 
postulated  that  the  reduced  gradient  may  be  the  result  of  wave- 
induced  motions  and/or  spume  droplet  production  (e.g.,  see  the 
discussionby  Wm  [1990]  mdDeLeeuw  [1990]).  These  spume 
drops  are  generated  when  spray  is  directly  tom  off  the  wave 
crests  in  high-wind  conditions  and  are  addressed  in  section  6. 

Finally,  in  the  Q.USE  setup  a  droplet  ejected  at  the  tunnel 
entrance  and  kept  aloft  for  the  entire  length  of  the  test  section 
(30  m)  has  approximately  4  s  (at  a  nominal  wind  speed  of  7.5 
m  s'')  to  interact  with  the  turbulent  fields  of  temperature  and 
humidity.  The  work  of  Andreas  [1990]  has  shown  that 
droplets  having  an  initial  radius  smaller  than  20  pm  are  able 
to  exchange  68%  of  their  mass  with  the  environment  in  4  s. 
We  can  see  the  result  of  this  exchange  in  Figure  10,  where  the 
contribution  of  the  droplets  to  the  moisture  field  increases 
significantly  with  fetch.  However,  this  also  means  that  much 
less  than  half  of  the  droplet  volume  (found  by  integrating  under 


the  curves  in  Figure  6)  has  a  chance  of  contributing  to  the 
moisture  field. 

6.  Marine  Boundary  Layer  Simulation 

The  above  discussion  is  particularly  relevant  to  open  ocean 
research  because  we  have  evidence  that  the  source  of  droplets 
generated  at  IMST  is  representative  of  high-wind  conditions 
over  the  open  ocean.  This  is  based  on  the  comparison  between 
tile  QJJSE  source  function  and  the  estimates  of  the  over-ocean 
function  given  by  Andreas  [1992]  shown  in  Figure  1 1,  which 
attempts  to  include  the  additional  source  of  droplets  arising 
from  spume  drop  production.  Figure  11  shows  that  the  two 
functions  are  comparable  for  wind  speeds  between  15  and  18 
m  s"'.  The  review  hy  Andreas  etal  [1995]  placed  this  function 
on  the  high  end  of  the  most  probable  estimates  of  spray  droplet 
production.  Therefore  this  function  is  expected  to  maximize 
the  influence  of  spray  droplets  on  the  near-surface  energy 
budget,  althougli  even  this  assumption  is  far  from  certain  due 
to  our  lack  of  knowledge  concerning  the  generation  of  the 
largest  (spume)  drops. 

In  order  to  examine  the  influence  of  the  droplets  under  high- 
NVind  conditions  in  the  MBL,  we  conducted  two  separate 
simulations  of  a  fully  developed  atmospheric  surface  layer 
using  the  Andreas  [1992]  source  function.  In  both  simulations 
we  initialize  Gwaihir  with  a  constant  flux  layer  using  a  lO-m 
wind  speed  J/jq  of  18  m  s'',  an  air  temperature  of  20®C,  a 
water  temperature  of  22 '"C,  and  a  relative  humidity  of  80%. 
Periodic  boundary  conditions  are  used  to  allow  the  droplets  to 
interact  with  the  turbulent  fields  as  long  as  they  are  airborne. 
The  velocity  is  fixed  at  the  upper  boundary,  while  the  tempera¬ 
ture  and  humidity  values  are  allowed  to  adjust  so  that  the  flux 
remains  constant  across  the  upper  boundary.  These  conditions 
are  specifically  chosen  for  comparison  with  the  analytical 
model  results  given  by  Andreas  etal.  [1995], 


Figure  11.  A  comparison  of  the  Gwaihir  source  function  with 
the  parameterization  of  the  over  ocean  source  function  given  by 
Andreas  [1992]. 
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Figure  12*  The  model-derived  sensible  heat  fluxes  computed 
in  the  absence  and  presence  of  droplets  released  at  two 
difference  locations.  The  conditions  during  the  run  were  a 
wind  speed  of  18  m  s*'  and  air  temperature  of  20°C  at  10  m,  a 
water  temperature  of  22  ®C,  and  a  relative  humidity  of  80%, 


The  current  version  of  Gwaihir  does  not  attempt  to  directiy 
model  the  influence  of  waves.  However,  we  believe  we  can 
place  some  quantitative  limits  on  the  role  they  play  in  droplet 
dispersion  and  evaporative  processes  through  the  height  at 
which  we  release  the  droplets.  In  the  first  simulation  we 
release  the  droplets  at  their  ejection  height,  while  in  the  second 
simulation  we  release  the  droplets  at  the  significant  wave 
height  determined  by 

0.015  C/fo  (44) 

The  latter  is  intended  to  simulate  a  situation  where  the  spume 
droplets  are  sheared  off  the  wave  crests  and  immediately  find 
themselves  at  a  considerable  distance  from  the  mean  surface. 
Andreas  [1992]  has  argued  that  is  an  appropriate 

timescale  in  his  model,  both  because  of  the  spume  droplet 
effect  and  because  it  serves  as  a  means  to  parameterize  the 
increase  of  turbulence  intensity  with  wind  speed.  However, 
because  turbulent  dispersion  is  already  included  in  our  model, 
we  expect  this  simulation  to  give  a  best  case  scenario  for 
droplet  dispersion.  This  simulation  should  therefore  approxi¬ 
mate  the  upper  bound  on  the  possible  droplet  influence  under 
these  conditions.  The  release  of  the  droplets  from  their  ejection 
height  is  far  less  favorable  for  droplet  dispersion  and  is 
expected  to  place  a  lower  bound  on  the  interactive  processes 
for  this  particular  source  function. 

The  droplet's  influence  on  the  latent  and  sensible  heat  fluxes 
for  the  two  runs  are  shown  in  Figures  12  and  13,  where  the 
fluxes  are  calculated  from  (9)  and  (10).  The  change  to  the 
mean  humidity  and  temperature  profiles  are  also  shown  in 
Figures  14  and  15  for  comparison  with  the  LBL  results. 


Figure  13.  The  model-derived  latent  heat  fluxes  computed  in 
the  absence  and  presaice  of  droplets  released  at  two  difference 
locations.  The  conditions  during  the  run  are  as  in  Figure  12. 

Figures  12  through  15  show  that  the  influaice  of  droplet 
evaporation  on  the  mean  profiles  and  surface  heat  fluxes  is 
strongly  affected  by  the  height  at  which  the  droplets  are 
released. 

The  droplets  released  at  their  ejection  height  have  a  minimal 
impact  on  the  latent  and  sensible  heat  fluxes  and  produce 
changes  to  the  mean  jjroflles  that  differ  only  slightly  from  the 
laboratory  results.  However,  the  droplets  released  at  the  wave 
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Figure  14,  Profiles  of  the  change  in  the  model-derived 
temperature  profiles  due  to  the  presence  of  droplets  in  our 
marine  boundary  layer  (MBL)  simulation.  The  conditions 
during  the  run  were  a  relative  humidity  of  80%  and  a  10-m 
wind  speed  of  18  m  s‘\ 
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Figure  15.  Profiles  of  the  change  in  the  model-derived  specific 
humidity  profiles  due  to  the  presence  of  droplets  in  our  MBL 
simulation.  The  conditions  during  the  run  were  a  relative 
humidity  of  80%  and  a  10-m  wind  speed  of  18  m  s*^ 


height  (4.86  m  in  this  simulation)  significantly  increase  an 
already  substantial  latent  heat  flux  by  20%  at  the  top  of  model 
domain.  Tlie  relatively  weaker  sensible  heat  flux  is  decreased 
by  more  than  1(K)%.  These  finding  agree  very  well  with  the 
results  of  Andreas  et  at.  [1995],  which  found  that  the  total 
droplet  generated  heat  flux  was  roughly  30%  of  the  interfacial 
(bulk)  flux. 

An  interesting  feature  of  these  profiles  is  the  apparent 
asymmetry  of  the  droplet  feedback  effects  described  in  part  1 
and  by  Andreas  et  al  [1995].  These  effects  act  to  reduce  the 
latent  heat  flux  near  the  surface  and  increase  it  above.  Fairall 
et  al  [1995]  hypothesized  that  this  feedback  could  be  included 
in  the  model  of  Andreas  [1992]  by  assuming  that  only  a 
fraction  of  the  total  droplet  gen^ated  heat  fluxes  appears  above 
the  droplet  evaporation  layer.  The  asymmetry  shown  in 
Figures  12  and  13  suggests  that  this  fraction  is  greater  than  the 
50%  used  by  Fairall  et  al  [1995]. 

7.  Conclusions 

In  this  paper  we  have  shown  that  the  Lagrangian  model 
simulations  given  in  part  1  have  been  even  further  improved  by 
simply  in^roving  the  simulation  of  the  turbulent  fields  through 
which  the  droplets  are  dispersed.  In  the  field,  where  it  is 
difficult  to  measure  even  mean  profiles  near  the  surface,  we  can 
use  the  k-e  model  to  provide  the  required  Lagrangian  model 
parameters  using  measurements  from  the  surface  and  at  some 
reference  height  to  initialize  the  Eulerian  model.  The  droplet 
model  can  then  be  used  to  examine  the  influence  of  the  droplet 
evaporation  and  sensible  heat  release  on  the  surface  energy 
budget  using  a  variety  of  wind-dependent  source  functions  once 
reliable  field-based  source  function  estimates  are  available. 


Alternatively,  Gwaihir  could  be  used  with  the  few  near¬ 
surface  profile  measurements  that  are  available  to  estimate  the 
source  function.  This  would  be  accomplished  by  adjusting  the 
model  source  function  until  the  modeled  profiles  match  the 
measurements.  The  major  problem  with  this  approach  is  our 
lack  of  knowledge  of  how  spume  drops  are  produced  over  the 
ocean.  This  leads  to  a  great  deal  of  uncertainty  in  how  to 
parameterize  their  production  in  numerical  models.  However, 
current  research  efforts  focusing  on  spume  drop  production 
should  improve  these  parameterizations  in  the  near  future. 

Our  results  indicate  that  an  increase  in  turbulence  intensity 
due  to  high  winds  does  not  significantly  increase  the  effect  that 
evaporating  jet  droplets  have  on  the  scalar  fields  of  temperature 
and  humidity.  The  principal  reason  for  this  is  that  the  turbu¬ 
lence  is  still  not  strong  enough  at  the  droplet’ s  ejection  height 
to  overcome  the  fall  velocity  of  the  largest  droplets.  However, 
the  potential  for  substantial  modification  of  the  surface  heat 
fluxes  exists  if  the  presence  of  waves  acts  to  eject  the  droplets 
higher,  permitting  the  larger  droplets  to  remain  aloft  for  longer 
poiods  of  time.  Therefore  we  need  to  gain  a  better  understand¬ 
ing  of  how  the  droplets  are  generated  at  high  wind  speeds  (i.e., 
as  spume  and/or  jet  drops),  as  well  as  how  the  waves  affect  the 
dispersion  of  these  droplets  once  airborne. 

An  initial  attempt  to  study  the  waves  effect  has  been  reported 
by  Andreas  et  al.  [1995]  using  modifications  to  the  model 
described  in  part  1 .  Although  the  relevance  of  these  results  as 
limited  by  the  simplified  flow  field  used  to  model  the  wave- 
induced  velocity  perturbations,  these  results  clearly  indicated 
that  the  presence  of  waves  significantly  increases  the  amount  of 
vertical  dispersion.  Therefore  the  logical  next  step  in  the 
development  of  Gwaihir  is  to  adapt  it  more  fully  to  over-ocean 
conditions  byintroducing  a  wavy  lower  surface  and  using  the 
model  to  simulate  the  flow  field.  We  expect  to  produce  a  more 
accurate  simulation  of  the  flow  over  waves  by  using  an 
approach  similar  to  the  modeling  work  of  Gent  and  Taylor 
[1976]. 

Gwaihir  will  also  have  to  determine  a  way  to  simulate  the 
release  of  spume  droplets  if  we  wish  to  accurately  model  the 
effect  of  spray  droplets  under  hi^-wind  conditions.  This  might 
involve  testing  proposed  source  functions  of  spume  droplets 
and  releasing  these  droplets  only  at  the  wave  crests.  Of  course, 
the  inclusion  of  these  droplets  leads  to  a  numb^  of  yet  unre¬ 
solved  questions  about  their  generation,  such  as  the  appearance 
of  flow  separation  over  breaking  waves.  However,  we  will 
leave  those  topics  for  future  discussion. 

Finally,  the  effects  introduced  by  modeling  droplets  com¬ 
posed  of  seawater,  rather  than  fresh  water,  would  have  to  be 
included.  This  would  primarily  involve  changes  to  the  equa¬ 
tions  governing  evaporation  rate  and  determination  of  how  to 
deal  with  the  residual  sea-salt  droplet  that  remains  behind  after 
evaporation.  Such  droplets  could  have  long  residence  times 
and  are  of  importance  in  processes  involving  cloud  physics  and 
atmospheric  optics, 
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1,  INTRODUCTION 

Ocean  scientists  have  been  speculating  for  about  50  years 
that  sea  spray  droplets  can  enhance  the  air-sea  fluxes  of  heat 
and  moisture.  The  Iterature  in  this  field,  however,  has  yield^ 
no  consensus  one  way  or  the  other.  One  of  the  major 
problems  involves  the  uncertainty  surrounding  the  source 
function  for  these  droplets.  The  number  of  droplets  produced 
at  the  sea  surface  is  known  to  be  a  function  of  both  droplet 
size  and  environmental  conditions.  However,  the  various 
source  functions  proposed  in  the  Iterature  spans  several  order 
of  magnitude..  Once  the  spray  droplets  are  airborne,  the 
principal  concern  of  modelers  is  to  determine  a  reaistic  means 
to  simulate  the  interaction  between  the  evaporating  droplets 
and  the  temperature  and  moisture  fields. 

Droplets  at  the  small  end  of  the  size  spectrum  remain 
suspended  indefinitely  and  thus  can  exchange  all  the  heat 
and  water  vapor  they  have  available.  But  because  of  the 
number  produced  and  their  size,  these  small  droplets  do  not 
contribute  much  to  the  air-sea  heat  and  moisture  fluxes.  Large 
droplets,  on  the  other  hand,  carry  a  lot  of  heat  and  moisture 
but  fall  back  into  the  sea  so  quickly  that  they  do  not  contribute 
much  to  the  surface  fluxes  either.  The  mid-range  droplets 
-  -  those  with  radii  at  formation  between  1 0  and  200  pm  (mostly 
spume  droplets)  -  thus  contribute  most  to  the  air-sea  fluxes 
because  they  are  small  enough  to  remain  suspended  for  a 
few  seconds  and  also  carry  a  moderate  amount  of  heat  and 
moisture. 

IndK/iduaHy.  we  have  developed  two  distinct  models  for 
the  role  of  sea  spray  in  air-sea  heat  and  moisture  exchange. 
Details  about  the  models  can  be  found  in  Andreas  (1990, 
1992),  Edson  and  Fairall  (1994),  and  Edson  et  al.  (1996). 
Therefore,  we  only  provide  a  brief  description  of  our  models 
in  the  following  section  and  then  focus  on  two  test 
comparisons.  These  two  comparisons  are  chosen  to  simulate 
conditions  where  we  might  expect  significant  spume  droplet 
production:  a  moderate-wind-speed  case  (10  m/s)  over  the 
tropical  Pacific,  and  a  high-wind-speed  case  (18  m/s)  at 
mid-latitudes. 

2.  THE  MODELS 

Andreas’s  (1992)  model  has  three  components.  First, 
it  predicts  that  rate  at  which  sea  spray  droplets  with  initial  radii 
between  2  and  500  pm  are  produced  at  the  sea  surface  as 
a  function  of  wrind  speed.  This  spray  generation  function  is 
derived  from  Millet's  (1 987)  generation  function  but  also  has 
the  first  realistic  prediction  for  spume  production.  Spume 
droplets  are  those  tom  directly  off  the  wave  crests  by  the  wrind , 
are  typically  20  pm  in  radius  and  larger,  and  contribute  most 
to  the  spray  sensible  and  latent  heat  fluxes  (And  reas  1 992) . 


The  second  component  of  Andreas’s  (1 992)  spray  model 
is  a  complete  mictophysical  model  that  computes,  for  dropists 
of  arbitrary  size,  time  scales  that  quantify  how  rapidly  individu^ 
droplets  exchange  sensible  and  latent -heat  with  their 
environment  (Andreas  1 990).  By  comparing  these  time  scales 
with  an  estimate  of  a  droplet’s  residence  time  above  the  sea 
surface,  the  model  can  estimate  how  much  of  a  droplefs 
available  heat  and  water  it  can  exchange  before  falling  back 
into  the  sea.  The  third  component  of  Andreas's  spray  model 
is  thus  an  estimate  of  this  residence  time,  parameterized  as 
the  quotient  of  the  significant  wave  ampitude  and  the  droplefs 
settfing  speed  in  still  air. 

Edson’s  model  consists  of  two  interactive  components 
(Edson  etal.,  1996).  The  first  component  is  an  Eulerian  k-e 
model  of  turbulent  flows  that  closes  the  system  of  equations 
using  two  prognostic  equations  forthe  turbulent  kinetic  energy, 
k,  and  its  rate  of  dissipation,  e .  The  second  component  is 
a  Lagrangian  model  of  evaporating  sea  spray  (Edson  and 
Fairall,  1994).  The  Lagrangian  model  keeps  track  of  the 
amount  of  sensible  heat  and  moisture  absorbed  or  released 
in  each  cell  of  the  Eulerian  domain  during  each  step  of  a 
droplet’s  trajectory.  This  determines  the  source/sink  ternis 
forthe  prognostic  equations  of  temperature  and  moisture  in 
the  Eulerian  model.  The  two  models  are  then  run  iteratively 
until  a  steady  state  is  reached. 

3.  RESULTS 

In  this  initial  comparison,  Edson’s  model  attempts  to 
include  the  presence  of  waves  by  simply  releasing  the  droplets 
at  the  same  significant  wave  amplitude  used  by  Andreas 
(1992).  The  model  also  imposes  constant  flux  conditions  at 
the  top  of  the  model  domain,  periodic  lateral  boundary 
conditions,  and  constant  surface  roughness  at  the  lower 
boundary.  The  same  source  function  is  used  in  both  models. 

Edson’s  model  is  able  to  predict  the  sensible  and  latent 
heat  flux  profiles  in  the  absence  (H,  and  H  J  and  presence 
(F,  and  F^)  of  spray  droplets.  The  profiles  of  H,  ar^  H  l  are 
shown  by  the  sold  ines  in  Figs.  1  and  2,  while  profiles  of  F, 
and  Fl  are  shown  by  the  broken  ines  forthe  two  simulations. 
Andreas’s  model  computes  H.  and  H^  and  the  nominal  spray 
sensible  (QJ  and  latent  (QO  heat  fkjxes.  The  two  models  can 
be  compared  using 

F,  =  H,  +  PQ,-aQ,-YQi  (1) 


F,=  H,-aQ,  (2) 


where  a,  P,  andy  are  small,  non-negative  numbers.  The  y 
term  accounts  for  possible  feedback  effects.  We  should  also 
point  out  that  Andreas  (1992)  computes  as  negative 
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because  droplets  must  consume  latent  heat  to  evaporate. 
Andreas  (1 992)  Implicitly  assumed  a  =  p  =  1,Y  =  0.  Fairall  et 
al.  (1994)  tried  a  =  P  =  0.5,  Y  =  Oln  a  model  examining  spray 
droplet  effects  in  tropical  storms.  We  examine  these  particuiar 
constants  In  the  comparisons  summarized  by  Table  1. 

4.  CONCLUSIONS 

Our  two  models  compare  reasonably  well  in  both  cases. 
In  particular,  both  predict  that,  for  the  high  wind  speed  case, 
the  spray  latent  heat  flux,  especially,  will  be  at  least  1 0%  of 
the  normal  turbulent  or  interfadal  latent  heat  flux.  In  other 
words,  sea  spray  can  significantly  enhance  the  air-sea  heat 
and  moisture  fluxes  in  high  winds.  Our  modeEng  also  suggests 
how  to  partition  the  total  air-sea  sensible  and  latent  heat  fluxes 
into  contributions  from  turbulent  and  spray  fluxes.  We  are 
now  working  on  the  inclusion  of  a  wavy  lower  boundary  in 
the  Eulerian  code.  We  would  then  release  the  droplets  at 
their  (jet  droplet)  ejection  height  along  the  wave  to  obtain  a 
more  realistic  simulation  of  the  influence  of  the  waves.  If 
succesdul,  we  beieve  that  we  will  be  able  to  devebp  a  simple 
parameterization  of  the  spray  heat  fluxes  for  use  In  large-scale 
models. 

This  work  is  supported  by  the  Office  of  Naval  Research. 
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Fig.l  Profiles  of  the  sensible  heat  flux  computed  In  the  absence 
(solid  line)  and  presence  (broken  line)  of  droplets. 


Fig,2  Profiles  of  the  latent  heat  flux  computed  in  the  absence 
(solid  line)  and  presence  (broken  line)  of  droplets. 


Table  1 .  Comparison  of  our  model  results  under  enviromental  conditions  commonly  experienced  during  autumn  in 
Mid-Latitudes  and  during  westerly  wind  bursts  in  the  Tropical  Pradfic. 


Flux  Variable 

Case  1 :  “Mid-Latitude" 

Case  2:  “Tropical  Padfic” 

Andreas 

Edson 
(a,P=1,  Y=0) 

Edson 

(a,3=.5,  Y=0) 

Andreas 

Edson 

(a.p=1,Y=0) 

Edson 

{a.p=.5.  Y=0) 

U,o  =18  m/s;  T*  =20°C;  T...  =22°C;  RH  =80% 

U,o  =10  m/s;  T*  =27»C;  =29"C;  RH  =80% 

H,  (W/m^) 

41.7 

44.9 

44.9 

26.5 

28.3 

28.3 

Hl  (W/m^) 

243.4 

250.1 

250.1 

224.3 

235.5 

235,5 

Q,  (W/m^) 

6.9 

1.6 

3.2 

0.9 

0 

0.1 

Q,  (W/m^*) 

-62.4 

-65.9 

-131.8 

-4.4 

-2.6 

-5.2 
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1.  INTRODUCTION 


(2)  must  be  modified  to  take  into  account  the  energy  into  the 
ocean,  F(0) ,  such  that 


Investigations  of  atmospheric  turbulence  over  the  world’s 
oceans  have  shown  that  the  interaction  of  wind  with  surface 
waves  results  in  flow  characteristics  that  differ  substantially 
from  even  a  horizontally  homogeneous  surface  layer  over  land. 
As  a  result,  marine  meteorobgists  and  physical 
oceanographers  often  divide  the  boundary  layer  close  to  the 
surface  into  the  surface  layer  where  Monin-Obukhov  (MO) 
similarity  holds,  and  a  wave  boundary  layer  (WBL)  where 
additional  scairng  parameters  are  required  for  similarity.  Even 
though  the  search  for  these  scaling,  parameters  and 
hypotheses  for  their  use  has  been  going  on  for  many  years 
(e.g..  Miles,  1957;  Hare  et.  al,  1997),  we  are  still  a  long  way 
from  a  consensus  in  the  sdentific  community. 

In  this  paper  we  present  some  of  our  on-going 
investigations  of  turbulence  and  wave-induced  flow  in  the 
marine  surface  layer.  These  investigations  rely  on  a  set  of  data 
collected  from  the  R/P  FLIP  during  the  Marine  Boundary  Layers 
(MBL)  experiment  sponsored  by  the  Office  of  Naval  Research. 
Spedfically,  this  paper  examines  the  energy  flux  and  Its 
relatbnship  to  its  rate  of  dissipatbn  within  the  marine  boundary 
layer.  The  flux  of  kinetb  energy,  F(z),  into  a  layer  of  air  over 
a  horizontally  homogeneous  surface  Is  given  by 


F{2)  -  uV'(y(z)  +  Vlfz)  +  +  ^w'p' 

P 


(1) 


h 


jedz  = 
0 


h 


-lF{h)-F{0)]--^[wXdz 

^0 


(3) 


Therefore,  less  volume  averaged  dissipation  is  required  to 
balance  the  same  energy  flux  into  the  layer  as  long  as  there 
is  a  net  flux  into  the  ocean.  In  fact,  the  difference  between 
the  integrated  dissipation  rate  and  the  total  energy  flux  at  the 
top  of  this  layer  can  be  used  to  roughly  estimate  the  energy 
flux  into  the  ocean. 


2.  SIMILARITY  THEORY 


In  a  stationary,  horizontally  homogenous,  constant  stress 
layer,  the  vertical  derivative  of  the  energy  flux  takes  the  form 
of  the  familiar  TKE  budget  equation 


dZ  dZ  dZ  p  dZ 


MO  scaing  provbes  us  with  a  dimensionless  form  of  the  TKE 
budget  given  by 


where  the  first  two  terms  on  the  right-hand-side  represents 
the  flux  of  mean  flow  kinetic  energy,  and  the  second  two 
represent  the  rate  of  diffusion  of  kinetic  energy.  Over  land, 
we  often  assume  that  the  energy  flux  through  the  ground  is 
negigible,  such  that  the  flux  entering  the  layer  at  height  h  can 
be  related  to  the  total  rate  of  dissipation  within  the  layer  by 

h  h _ 

fedz  =  -F{h)  *  -^1  wXdz  (2) 

0  ''^0 


where  the  second  term  on  the  right-hand-side  accounts  for 
the  generation  of  kinetic  energy  due  to  any  buoyancy  flux. 
For  neutral  conditions,  this  expression  states  that  the  flux  of 
kinetic  energy  into  a  layer  is  balanced  by  the  total  rate  of 
dissipation  within  that  layer. 

Over  the  ocean,  we  expect  the  surface  energy  flux,  whbh 
drives  waves  and  currents,  to  be  non-negligible.  Expression 
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where/. istheMOtength,  k  isthevonKarman constant, and  u, 
is  the  friction  velocity.  This  expressbn  is  often  used  to  estimate 
the  momentum  over  the  ocean  from  estimates  of  the  dissipation 
rates  and  a  parameterization  of  from 

puf  =  p[eKZ/(t)^(2/L)p  (6) 


This  type  of  parameterization  should  be  valid  as  long  as 
the  measurements  are  made  above  the  WBL  and  the  constant 
flux  assumption  holds.  However,  an  additional  forcing 
mechanism  is  introduced  due  to  the  presence  of  su  rface  waves 
as  we  approach  the  surface.  In  this  region  (i.e.,  the  WBL), 
MO  similarity  is  expected  to  break  down.  For  example,  in  the 
WBL  the  velocity  can  be  decomposed  Into  mean,  turbulent, 
and  wave-induced  components 

u(t)  =  U-Fu'{()^Q(t}  (7) 
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where  the  wave-induced  component  is  defined  using  an 
extensbn  of  Reynolds  averaging  known  as  phase  averaging 
(e.g,,  Finnegan  et  al. ,  1984). 

This  type  of  averaging  can  be  used  to  rewrite  the  TKE 
energy  budget  to  include  the  wave-induced  components.  For 
example,  the  shear  production,  P,  would  include  terms 
representing  the  energy  production  due  to  the  interaction 
between  the  wave-indued  flux  and  the  mean  shear 

P=-uV[l  -  Vw'[\ 

'■  dz  dz'  ' 


The  bracketed  terms  are  not  expected  to  obey  M-0  scaling, 
which  leads  to  a  breakdown  of  similarity  relationships  like  (5) 
within  the  WBL 

3.  RESULTS  AND  CONCLUSIONS 

Even  if  we  are  able  to  separate  the  wave-induced 
fluctuations  from  the  turbulence,  a  scale  analysis  of  even  a 
simplified  form  of  kinetic  energy  equation  is  a  difficult  task. 
However,  we  can  begin  to  make  some  inferences  about  the 
behavior  of  the  combined  fluctuations  and  their  related  fluxes 
(i.e.,  the  total  flux)  using  the  data  collected  during  the  MBL 
experiment.  The  data  collected  from  the  FLIP  mast  allows 
us  to  compute  all  the  terms  of  the  total  kinetic  energy  flux  at 
2-3  levels,  and  the  dissipation  rate  of  TKE  at  4  levels. 

We  begin  our  discussion  with  a  comparison  of  our  e 
measurements  with  their  MO  similarity  prdiction 

(9) 


where  we  assume  a  balance  between  production  and 
dissipation  for  4)^,  i  e.,  we  have  not  attempted  to  parameterize 
the  transport  terms.  We  expect  (9)  to  be  greater  than  our 
measured  values  as  bng  as  there  is  a  net  flux  Into  the  ocean. 
Evidence  for  this  is  shown  in  Fig.  1 ,  which  shows  that  the 
measured  dissipation  profiles  are  less  than  the  M-0  prediction. 
These  profiles  were  averaged  over  a  period  of  increasing  winds 
and  growing  seas.  The  “dissipation  defidt”  becomes  smaller 
as,  presumably,  we  move  out  of  the  WBL. 

We  attempt  to  recondle  the  difference  between  the 
predicted  and  measured  dissipation  rates  by  considering  the 
underlying  sea-state  in  Fig.  2.  In  this  figure,  we  have 
normalized  our  dissipation  estimates  by  (9)  and  then  bin- 
averaged  according  to  the  indicated  ranges  of  the  wave-age 
parameter,  c/u,;and  /cz,  where  k  is  the  wavenumber  of  the 
dominant  waves  defined  from  the  peak  of  the  wave  spectra. 

Our  dissipation  estimates  are  adequately  predicted  by 
MO  similarity  theory  for  mature  seas  above  /cz  «  1 .5 .  However, 
our  results  also  indicate  that  their  is  an  appredable  difference 
between  our  dissipation  measurements  and  their  MO 
predictions  for  all  sea-states  when  kz<^ .  The  dissipation 
deficit  is  greatest  over  the  youngest  sea  as  one  would  intuitively 
expect.  Therefore,  our  results  indicate  that  the  dissipation 
rate  does  not  obey  M-0  similarity  within  the  WBL,  and  that 
the  dissipation  deficit  appears  to  be  a  function  of  sea-state. 

We  are  now  inves^ating  whether  the  dissipation  deficit 


is  balanced  by  the  transport  of  energy  to  the  surface  and/or 
by  a  modification  of  the  energy  production  due  to  the  presence 
of  waves  (see  abstracts  by  Hare  et  al.,  Hristov  et  a!.,  and  Miller 
et  al.  in  this  volume).  Our  findings  have  important  impications 
foruseof  the  inertial  dissipation  and  buk  aerodynamic  methexis 
over  the  cx^ean,  particularly  from  buoys.  The  influence  of 
the  surface  waves  on  the  near  surface  turbulence  is  significant 
and  a  means  to  account  for  their  presence  should  be  included 
in  subgrid-scale  parameteiizations  in  LES  and  bwer  boundary 
conditions  in  numerical  models. 

This  work  was  supported  by  the  Office  of  Naval  Research. 
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Fig.  1 .  Measured  and  predicted  dissipation  profiles  averaged 
over  a  4  hour  period. 
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Fig.  2.  Bin-averaged  normalized  dissipation  profiles. 
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1  INTRODUCTION 

Recent  advancements  in  the  technology  of  the  experiment 
made  it  possible  to  collect  high  quality  field  data  sets,  which 
contain  better  information  about  the  coupling  between  the 
wind  and  the  waves  and  can  be  used  to  approach  in  a  new 
way  longstanding  problems  related  to  wind-waves  interac¬ 
tion.  All  this  stimulates  the  development  of  new  methods 
for  data  processing  and  analysis.  In  this  paper  we  present 
and  apply  such  a  new  method.  The  data  analyzed  here  are 
from  the  Marine  Boundary  Layer  Experiment,  which  took 
place  in  April-May  1995  50  kilometers  off  the  coast  of  Cali¬ 
fornia.  The  instrumentation  used  in  this  analysis  includes  a 
vertical  array  of  sonic  and  cup  anemometers  for  horizontal 
and  vertical  velocity  measurements,  and  a  wave  wire  situ¬ 
ated  directly  beneath  this  array  to  compute  the  necessary 
wave  statistics. 

The  air  flow  over  ocean  waves  is  conveniently  decom¬ 
posed  into  mean,  turbulent  and  wave  coherent  components 
u  =  u  -h  uturb  +  u.  The  separation  is  made  by  the  ori¬ 
gin  of  the  components,  but  also  more  importantly,  by  the 
fact  that  the  coupling  between  these  components  is  weak 
enough.  The  wave  coherent  component  is  the  one  responsi¬ 
ble  for  the  exchange  of  momentum  and  energy  between  the 
waves  and  air  flow.  That  is  why  its  identification  from  field 
experiment  data  is  of  primary  interest. 

To  approach  the  identification  problem  we  can  assume 
that  the  turbulent  component  in  the  air  flow  and  the  wave 
are  not  correlated.  Let  us  consider  monochromatic  waves 
with  a  period  To  and  the  air  flow  velocity  at  a  certain  height 
above  the  point  where  the  wave  height  is  measured.  Under 
the  assumptions  above,  averaging  the  air  flow  velocity  for  a 
number  of  periods  N, 


^To(i)  =  Su(t)  =  ^u(t  +  uTo)  j  -  tA,  (1) 

should  filter  out  the  turbulence  and  what  will  be  left  is  the 
phase  average  of  the  wave  coherent  component  [1,2].  If  the 
waves  are  monochromatic  such  a  technique  should  perform 
well.  However,  the  waves  in  the  open  ocean  are  spread  over 
an  interval  of  spectral  scales  and  the  attempt  to  use  (1)  with 
field  data,  leads  to  strong  attenuation  of  the  phase  average 
as  a  result  of  destructive  interference  [3].  This  motivated  the 
development  of  a  different  approach,  which  would  be  much 
less  sensitive  to  the  lack  of  coherence  in  nonmonochromatic 
wave  fields. 


2  PHASE  AVERAGING  VIA  THE  HILBERT  TRANS¬ 
FORM 

Just  to  introduce  the  idea,  let  us  recall  that  a  signal  which 
consists  of  two  Fourier  components  and  has 

instantaneous  amplitude  A{t)  and  phase  v?(t)  formed  by  the 
rule  For  a  signal  of  arbitrary 

spectrum  s{i),  the  idea  can  be  generalized  by  employing  the 
concept  of  the  analytic  signal  based  on  the  Hilbert  transform 
[4].  The  analytic  signal  ^(t)  is  a  complex-valued  function 
of  time  defined  as 


<!'(<)  =  +  (2) 

where  the  function  5H(f)  is  the  Hilbert  transform  of  5(t) 

SH{t)  =  Hs{t)  =  r  ^dr,  (3) 


and  V  indicates  that  the  integral  is  taken  in  the  sense  of 
the  Cauchy  principal  value.  From  (3),  the  Hilbert  transform 
SH{t)  of  5(t)  may  be  seen  as  the  convolution  of  the  functions 
5(t)  and  l/7rt.  For  the  frequencies  of  interest  in  physics 
w  >  0,  the  Fourier  transforms  of  s{t)  and  5H(f)  satisfy 
5^(a;)  =  — t5(u;);  i.e.,  ideally  snii)  may  be  obtained  from 
s{i)  by  a  filter  whose  amplitude  response  is  unity,  and  whose 
phase  response  is  7r/2  lag  at  all  frequencies  [4]. 

Let  T/(t)  be  the  wave  height  signal  over  a  time  interval  T 
and  let  u(t)  be  the  wind  velocity,  measured  above  the  point, 
where  the  wave  height  is  measured.  Let  ^(t)  = 
and  ^{t)  =  Arg^(t)  be  the  wave's  instantaneous  phase, 
as  defined  above.  Then,  to  select  and  average  the  wind 
velocities  at  the  moments  r  when  the  wave  phase  $(r)  = 
^(f),  we  can  use  the  transformation 


Vu{i)  = 


f^ujr)  ^r)  g($(r)-$(f))dr 
$'(r)  Smr)  -  m)  dr 


u. 


(4) 


For  monochromatic  signals  i/(t),  (4)  is  equivalent  to  (1). 
For  nonmonochromatic  exprerimental  signals,  S  and  P  show 
very  different  results  (Fig.  1). 


3  DISCUSSION 

We  can  use  Tu{t)  as  a  strictly  periodic  approximation  of 
the  wave  coherent  component  in  the  wind  u{t)  to  cal¬ 
culate  the  wave  coherent  momentum  and  energ;^  fluxes. 
Let  u{t)  =  {u(i),v{t)yw{t))  be  approximated  by  Vu{t)  = 
{uo  cos(wot  +  ^u),  Vo  cos{uot-\-  <l>v)i  tuo  cos(wot-f  )}.  and 
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Figure  li  A  wave  height  signal  vs,  time  (a),  its  syn^ronized 
average  5f/(t)  (b)  and  Hilbert-transform  phase  average  (c) 

(solid  lines).  In  (b)  and  (c)  the  dotted  lines  show  the  original 
signal.  The  period  of  the  most  energetic  mode  To  is  the  length 
of  the  subintervals  in  (b).  Due  to  destructive  interference,  the 
average  5T7(t)  of  10  periods  is  attenuated  and  distorted  (b).  In 
contrast.  Vrj{t)  on  (c)  shows  no  attenuation  or  distortions. 


the  wave  coherent  air  pressure  fluctuations  on  the  interface 
p{t)  by  Vp{t)  =  po  cos{uot  where  all  the  phases  are 

referenced  to  the  phase  of  =  r/ocoswot  and  uq  is 

an  appropriate  frequency.  Then  the  wave  coherent  stress 
(u{t)w{t))  and  the  energy  flux  (p(t)  dr\{i)ldi)  are 

Vw{t)^  =  0.5  UqWq  COs(<^u  4vi)  (5) 

(Vp{t)  aVr){t)ldt)  =  0.5  o^i]opo  sin((^p).  (6) 

The  results  (5)  and  (6)  make  obvious,  that  the  phase  dif¬ 
ferences  control  the  direction  of  the  momentum  and  energy 
fluxes  and  suggest  that  representing  the  structure  of  the  sur¬ 
face  layer  in  terms  of  amplitudes  and  j)hases  of  the  phase 
averaged  wave  coherent  fields  'Pu(t),  'Ppit)^  etc.,  would  be 
descriptive  and  adequate.  The  phase  shift  between  the  wave 
coherent  components  in  the  air  flow  and  the  wave  itself  in¬ 
dicates  a  separation  of  the  air  flow  [5,  6,  7,  8],  which  now 
can  be  detected  and  measured  by  using  the  transformation 
V  from  open  ocean  experiment  data.  Figure  2  shows  two 
results  obtained  with  the  phase  averaging  technique  we  just 
described. 

4  CONCLUSION 

In  conclusion,  we  presented  a  robust  approach  allowing  to 
separate  the  turbulent  and  wave  coherent  components  in  the 
air  flow  over  surface  waves.  The  technique  makes  it  possible 
to  identify  flow  separation  and  to  estimate  the  exchange  of 
momentum  and  energy  between  the  air  flow  and  the  waves 
from  field  experiments.  The  method  also  allows  to  detect 
inverted  momentum  and  energy  fluxes  for  decaying  seas  —  a 
phenomenon  predicted  by  models,  but  not  observed  in  open 
ocean  experiments. 


Figure  2:  The  left  plot  represenU  the  profiles  of  the  wave- 
coherent  stress  as  a  fraction  of  the  total  stress  -  (uti))/uj.  Each 
profile  coresponds  to  a  different  wave  age  Cp/UiOt  {Op  is  the 
waves'  phase  speed  and  C7io  is  the  mean  wind  speed  at  10  meters) 
as  shown  on  the  legend.  The  sign  of  the  stress  is  chosen  to 
be  positive  for  momentum  flux  downwards.  The  profiles  confirm 
that  for  developing  seas  most  of  the  wave  coherent  stress  occurs 
dose  to  the  surface  and  is  positive.  For  large  wave  ages  the  wave 
induced  stress  represents  a  considerable  fraction  of  the  total  stress 
and  is  negative  (upwards).  The  plot  on  the  right  shows  the  phase 
offset  of  the  pressure  fluctuations  vs.  wave  age  Cp/UiQ.  For 
low  wave  ages  the  pressure  is  lagging  the  wave  (the  duster  of 
positive  pressure  phases),  so  the  waves  grow.  For  larger  wave 
ages  (waves  out-running  the  wind),  the  pressure  is  leading  the 
waves  (higher  pressure  before  the  wave),  which  makes  the  waves 
to  decay.  The  sensor  measuring  the  pressure  fluctuations  was 
positioned  12  meters  above  the  surface. 
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6.  Continuing  work  on  the  simulation  of  the  flow  of  air  over  ocean  waves. 

The  Langrangian  and  coupled  models  have  been  validated  in  the  large  air-sea  interaction 
simulation  tunnel  at  IMST,  Luminy,  France.  As  such,  the  model  has  been  shown  to  perform  well 
in  simulations  of  a  developing  boundary  layer  over  a  roughened  surface.  Since  that  time,  the  Pis 
efforts  have  focused  on  inclusion  of  a  wavy  lower  boundary  layer  in  the  Eulerian  code.  This 
section  summarizes  these  efforts. 

Spray  droplet  modeling  3.  Simulating  droplet  dispersion  over  ocean  waves 

James  B.  Edson 

Woods  Hole  Oceanographic  Institution 


1.  Introduction 

Spray  droplets  are  produced  over  the  ocean  by  bubbles  bursting  at  the  surface.  These 
bubbles  are  usually  produced  by  the  air  entrained  by  breaking  waves.  These  bubble-generated 
spray  droplets  generally  fall  under  two  categories:  film  droplets  that  arise  from  the  breakup  of 
the  leading  side  of  the  bubble  (i.e.,  the  bubble  film)  as  it  penetrates  the  surfaee,  and  jet  droplets 
that  result  from  the  breakup  of  the  column  of  water  (i.e.,  the  jet)  ejected  from  the  collapsing 
bubble  cavity.  This  production  mechanism  also  accounts  for  the  generation  of  spray  droplets 
produced  by  rain  striking  the  surface.  Sea  spray  known  as  spume  droplets  are  also  produced 
through  direct  shearing  of  the  wave  erests  under  high  wind  conditions.  Therefore,  except  in 
instances  of  rain  under  light  winds,  these  spray  droplets  are  ejected  into  the  wave  boundary  layer 
(WBL)  where  wave-induced  motions  are  a  significant  component  of  the  velocity  field.  As  such. 


24 


the  dispersion  of  these  droplets  is  expected  to  be  influenced  by  the  presence  of  waves. 


2.  Model  Summary 


In  developing  equations  designed  to  study  flows  where  the  mean  departure  from 
hydrostatic  equilibrium  can  be  nonzero  (e.g.,  around  a  building),  Sini  [1986]  and  Sini  and 
Dekeyser  [1987]  decomposed  this  departure  from  hydrostatic  equilibrium  into  mean  and 
fluctuating  parts.  This  approach  is  also  well-suited  for  simulating  the  airflow  over  waves  as  the 
pressure  field  is  expected  to  be  out  of  local  hydrostatic  equilibrium.  Application  of  the 
Boussinesq  approximation  results  in  Reynolds  averaged  equations  for  the  mean  variables  given 
by 
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where  Einstein's  summation  notation  is  used;  the  overbar  represents  an  ensemble  average;  lower 
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case  letters  denote  the  turbulent  fluctuations;  v  is  the  kinematic  viscosity;  0^,  is  the  ambient 
virtual  potential  temperature;  0^  is  the  virtual  potential  temperature  of  the  reference  state  of  the 
fluid;  Q  is  the  total  specific  humidity;  g.  =  (0,0, -g)  where  g  is  the  gravitational  acceleration; 
is  the  density  of  air;  is  the  specific  heat  at  constant  pressure;  and  5^  and  represent  source 
terms  for  sensible  heat  and  moisture,  respectively.  These  source/sink  terms  simulate  the 
exchange  of  sensible  and  latent  heat  between  the  evaporating  droplets  and  the  scalar  fields  as 
they  are  dispersed  in  the  surface  layer. 

2.1.  Closure 


The  Reynolds-stress  tensor  is  modeled  using  a  slightly  modified  form  of  the  Boussinesq 
eddy  diffusivity  concept  to  allow  for  anisotropy  (see  Taylor,  1977  for  the  isotropic  forms) 
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where  a^,  a^,  anda^^  are  the  standard  deviations  of  the  longitudinal,  lateral,  and  vertical  velocity 
fluctuations,  respectively;  v^.  is  the  eddy  diffusivity;  and  e  is  the  turbulent  kinetic  energy  (TKE) 
defined  as 
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(9) 
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In  first-order  closure  models  the  eddy  diffusivity  in  a  neutral  surface  layer  is  parameterized  as 


(10) 


where  is  the  friction  velocity  andK  is  the  von  Karman  constant,  which  we  assign  the  value 
0.4.  In  our  e-e  closure  model,  the  eddy  diffusivity  is  parameterized  as 


(11) 


where  e  is  the  dissipation  rate  of  TKE  and  is  a  model  coefficient.  The  model  coefficient  can 
be  determined  using  the  near-surface  relationship  for  the  dissipation  rate  under  neutral  conditions 
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which  results  in 
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(13) 


Our  modifications  to  the  Boussinesq  concept  approximates  the  empirical  results  given  in 
Panofsky  and  Dutton  (1984)  such  that 
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Using  these  values  the  numerical  constant  becomes  =  1/36.  The  scalar  fluxes  are  then 
modeled  using 
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(15) 
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where  the  Prandtl  and  Schmidt  numbers  for  turbulent  diffusion  are  assigned  the  same  value, 

Pvj.  =  Scj  =0.95  [Hogstrom,  1988]. 

Closure  is  then  accomplished  through  prognostic  equations  for  both  the  TKE  and  its  rate 
of  dissipation 
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where  o^,  C^j,  o^,  and  are  additional  model  coefficients.  These  coefficients  are  derived 
using  the  expressions  in  Edson  et  al.  (1996)  that  provide  the  values  summarized  in  Table  1.  The 
production  terms  in  the  above  equations  are  also  modified  by  the  inclusion  of  terms  that  normally 
cancel  by  invoking  continuity.  For  example,  the  production  term  in  the  TKE  equation  can  be 
written  out  as 
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where  the  last  term  on  the  right  hand  side  is  normally  neglected  when  the  isotropic  form  of  the 
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Reynold  stress  tensor  is  used. 


Table  1.  Numerical  constants  used  in  present  model. 

Cel 

Ce2 

<Je 

K 

PVj 

Scj 

0.03 

1.44 

1.92 

1.0 

2.0 

0.4 

0.95 

0.95 

2.2  Model  ModiHcations  for  the  Undulating  Surface 

The  e-e  model  described  by  Eds  on  et  al.  (1996)  requires  a  number  of  modifications  to 
simulate  the  flow  over  ocean  waves.  An  obvious  difficulty  in  using  the  set  of  equations  given 
above  is  how  to  define  the  velocity,  temperature,  humidity  fields  when  z  is  beneath  the  wavy 
surface.  One  commonly  used  method  to  get  around  this  problem  involves  setting  up  a  new,  wave 
following  coordinate  system  ( X>  Z )  where 


X  =  X 

Z  =  Z  -  T)  e 


(20) 


where  q  is  the  instantaneous  wave  height  and  k  is  the  principle  wavenumber  of  the  undulating 
surface.  This  coordinate  system  has  the  desirable  characteristics  that  z  =  q  when  Z  =  0,  and 
z~  Z  for  Z  >  ^.  As  in  Taylor  (1977),  we  assume  that  the  flow  is  well  represented  by  logarithmic 
velocity  profiles  and  transform  Z  to  a  logarithmic  vertical  coordinate  C  as 


C  =  ln(Z) 


(21) 


The  set  of  governing  equations  can  be  rewritten  for  the  ( x ,  C)  coordinate  system  using  the 
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relationships 


^  +  j_  0^  az 

dx  dx  Z  dC  dx 


^  52 

dz  Z  dC  dz 


dx^  "  ax'  ^  Z  axac  a^  ^  z^  ac'  j  ^  z  ac  a;c2 


(24) 


a^(|)  ^  1  0^4)  ( azV 
az^  "  z^  ac"  i  az; 


^  1  0(1)  a"z 
^  z  ac  0^2 


(25) 


where  <J)  represents  the  appropriate  velocity,  temperature  or  humidity  variable  in  our  set  of 
equations. 

2.3  Initialization  and  Boundary  Conditions 

The  height  of  the  domain  is  chosen  such  that  =  3 Ik.  Periodic  boundary  conditions 
are  used  at  the  lateral  boundaries.  A  stokes  wave  is  chosen  for  the  lower  boundary 

H 

'n(X)  =  -^cos(A:x  -  (Of)  +  -^kcos(2ikx  -  (of)) 

2  8 

where  (o  is  the  angular  frequency  andH^  is  the  wave  height.  In  order  to  reach  steady-state 
conditions,  we  require  the  total  momentum  flux  to  be  constant  with  height  within  the  model 
domain  such  that 
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where  v  is  the  kinematic  viscosity  and  the  primes  and  tildas  denote  turbulent  and  wave  induced 
fluctuations,  respectively  [i.e.,  U(t)  -=U  +  u(0  =  U  +  u  \t)  +  u(t)].  The  last  term  on  the  right  hand 
side  of  (26)  represents  the  viscous  shear  stress  and  is  accounted  for  in  the  model  by  adding  v  to 
our  parameterization  of  Vj..  The  initial  wind  profiles  throughout  the  model  domain  are  defined 
in  a  frame  of  reference  moving  with  the  phase  speed  of  the  dominate  wave,  ,  such  that 


U(U)  -  +  —  ln| 

K 


C  +  C 


- 


(28) 


w(C,x) 


dT| 

dt 


-  {u.  *  i/.cx)) 


5r| 

dx 


(29) 


where  UJjO  =  A:t|(x)c^  is  the  wave  orbital  velocity,  and  is  the  surface  drift  current.  At  the 
upper  boundary  where  the  wind  speed  is  fixed  at 
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The  initial  simulations  rely  on  the  dispersion  relationship  for  deep  water  waves  to  provide 
the  relationships 
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The  wave  height,  phase  speed,  surface  roughness,  and  drift  current  are  parameterized  using 


=  0.015 
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where  t/jQ  =  t/(C  =  10)  and  P  is  the  Chamock  constant,  which  is  assigned  the  value  of  0.01 1. 
Therefore,  the  model  is  initialized  by  choosing  a  value  for  iteratively  solving  (31)  for  m,, 
and  computing  the  phase  speed  and  wave  height  from  U^q  =  +  ujK\n{\0l^^.  This  value 

of  the  friction  velocity  defines  the  constant  value  of  the  total  momentum  flux  throughout  the 
domain  as 


T  = 


(34) 


3.  Work  in  Progress 


The  remaining  component  of  the  model  that  needs  to  be  properly  parameterized  is  the 
energy  and  momentum  flux  boundary  conditions,  and  the  effects  of  the  wave-induced  energy  flux 
divergence.  The  momentum  flux  going  into  the  waves  (i.e.,  the  form  drag)  is  given  by 


-pMW^ 
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where  the  pressure-slope  correlation  represents  the  momentum  flux  supported  by  the  entire 
spectrum  of  waves.  The  present  model  can  only  resolve  the  form  drag  due  to  the  dominate 
wave.  Therefore,  the  question  arises  how  to  take  into  account  the  unresolved  component  of  the 
form  drag.  Our  model  simple  assumes  that  the  total  momentum  flux  minus  the  resolvable  form 
drag  must  be  balanced  by  the  shear  stress 
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The  mean  wind  at  the  first  grid  point  is  then  assigned  the  appropriate  value  to  provide  this 
balance. 

Another  difficulty  involves  the  non-zero  energy  flux  at  the  surface.  This  flux  is  mainly 
carried  by  the  wave  induced  pressure-vertical  velocity  correlation  at  the  surface 
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which  again  is  partially  resolved  by  our  model.  Therefore,  our  initial  simulations  assume  that  the 
total  energy  flux  minus  the  resolvable  flux  must  be  balanced  by  the  diffusive  flux 
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which  provides  the  value  of  the  kinetic  energy  at  the  lowest  grid  point.  The  total  energy  flux  is 
parameterized  using 

=  jyS(k)dk 

—  00 


where  yis  the  wave  input  term  (e.g..  Plant,  1982)  and5(^)  is  the  wave  spectum.  The  divergence 
of  this  flux  plus  the  wave  induced  energy  flux  appears  in  the  kinetic  energy  budget  as 

_d_ 
dz 

At  present,  only  the  divergence  of  the  turbulent  components  are  accounted  for  in  (17).  Therefore, 
proper  simulation  of  the  flow  over  waves  must  find  a  method  to  properly  account  for  these 
effects.  Several  methods  have  been  presented  in  the  literature  (e.g.,  Makin  and  Mastenbroek, 
1996;  Mastenbroek  et  al,  1996.  We  are  now  starting  to  compare  our  model  simulations  with 
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measurements  from  the  1995  Marine  Boundary  Layers  program  to  see  if  any  these  methods  or 

some  new  approach  can  provide  adequate  agreement  with  our  observations. 
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The  project  involved  an  investigation  of  the  effects  of  evaporating  sea-spray  under  high  wind  conditions. 
These  investigations  were  accomplished  using  a  interactive  Eulerian-Lagrangian  model  developed  by  the 
PL  The  interactive  model  was  validated  using  data  collected  during  the  1988  CLOSE  program.  This  work 
has  demonstrated  that  the  combined  model  accurately  simulates  the  turbulent  transport  of  evaporating 
droplets.  In  additions,  this  paper  advanced  the  state-of-the-art  in  droplet  research  by  demonstrating  that 
the  potential  for  substantial  modification  of  the  surface  energy  budget  exists  if  the  presence  of  waves  acts  to 
eject  the  droplets  higher  and/or  disperse  the  droplets  more  efficiently.  The  model  predicted  that  the  contri¬ 
bution  of  the  sea  spray  on  the  latent  heat  flux  is  at  least  10%  of  the  total  under  high  wind  speed  conditions. 
The  second  part  of  this  investigation  has  involved  modifications  to  the  Eulerian  portion  of  the  code  to 
include  a  wavy  lower  boundary.  The  validation  of  this  model  is  being  accomplished  through  comparisons 
with  an  open  ocean  data  set  collected  aboard  the  R/P  FLIP  during  the  1995  Marine  Boundary  Layers 
Experiment.  These  observations  are  now  being  used  to  improve  the  boundary  conditions  and  closure 
schemes  to  simulate  the  flow  over  waves  as  realistically  as  possible. 
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